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ABSTRACT 

Recent developments in the area of nonlinear control theory have shown how co- 
ordinate changes in the state and input spaces can be used with nonlinear feedback 
to transform certain nonlinear ordinary differential equations into equivalent linear 
equations. These feedback linearization techniques are applied to resolve two prob- 
lems arising in the control of spacecraft equipped with control moment gyroscopes 
(CMGs). The first application involves the computation of rate commands for the 
gimbals that rotate the individual gyroscopes to produce commanded torques on the 
spacecraft. The second application is to the long-term management of stored mo- 
mentum in the system of control moment gyroscopes using environmental torques 
acting on the vehicle. 

An approach to distributing control effort among a group of redundant actua- 
tors is described that uses feedback linearization techniques to parameterize sets of 
controls which influence a specified subsystem in a desired way. The approach is 
adapted for use in spacecraft control with double-gimballed gyroscopes to produce 
an algorithm that avoids problematic gimbal configurations by approximating sets 
of gimbal rates that drive CMG rotors into desirable configurations. The algorithm 
solves a normed approximation problem to match the desired rates using an affine 
parameterization of all gimbal rates that produce a specified acceleration of the 
spacecraft. 

The momentum management problem is stated as a trajectory optimization 
problem with a nonlinear dynamical constraint. Feedback linearization and collo- 
cation are used to transform this problem into an unconstrained nonlinear program. 
The procedure is developed in a general fashion for application to a broad class of 
linearizable systems. The approach to trajectory optimization is fast and robust to 
starting conditions in comparison to variational methods. A number of examples 
are presented showing applications to the proposed NASA space station. 

Thesis Advisor: Professor Derek Rowell 
Department of Mechanical Engineering 
Technical Supervisor: Edward Bergmann 
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Chapter 1 
Introduction 


The need to regulate the attitude of an orbiting spacecraft defines an important part 
of the overall design problem for a spacecraft control system. For example, the pro- 
posed NASA space station is required to provide stable inertial and Earth-pointing 
platforms for scientific instruments, and a platform capable of tracking objects. It 
is also necessary to prevent large attitude excursions during Space Shuttle docking 
and during re-boost maneuvers. The space station is being designed to orbit the 
Earth at a fixed attitude with respect to the local horizontal. In order to meet a 
wide variety of mission requirements, there are limits to the allowable deviations 
from the design attitude. These limits must be enforced even in the presence of 
uncertain environmental effects and on-board disturbances. Examples of environ- 
mental effects include aerodynamic and gravity torques, and on-board disturbances 
include crew motion, motion of the Mobile Manipulator, and the motion of the solar 
panels as they track the sun. As a consequence, it is necessary to develop algorithms 
which can utilize a set of available actuators to apply torques to the vehicle in order 
to maintain the attitude within the specified limits. 

There are a variety of actuators available to generate torques on a spacecraft for 
the purpose of attitude control. While propulsive actuators such as reaction control 
jets are probably the most common, there is another class of actuators which operate 
by transferring momentum to and from the spacecraft via a momentum ’’storage” 
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system. The two types of actuators which operate in this way are reaction wheels 
and control moment gyroscopes (CMGs). A reaction wheel consists of a flywheel 
whose axis of rotation is fixed with respect to the spacecraft. Torques are applied to 
the vehicle by speeding up or slowing down the spin rate of the flywheels. A CMG is 
constructed with a rotor that spins at a constant speed and is mounted in a gimbal 
system which can reorient the flywheel. CMGs can be used to apply torques to 
a spacecraft by rotating the angular momentum vectors of the individual rotors. 
Moving the rotors changes the net angular momentum of the group of CMGs. In 
the absence of external torques, momentum conservation implies a corresponding 
change in the vehicle momentum, or equivalently, the angular rate of the vehicle. 

There are a number of other devices which can be used to apply torques to 
orbiting spacecraft. Magnetic torquers produce torques on a vehicle through the 
interaction of the Earth’s magnetic field with currents in coils of wire. It is also 
possible to adjust panels attached to the spacecraft to alter the distribution of 
aerodynamic forces on the vehicle caused by the thin atmosphere encountered in 
low Earth orbit [1]. These last two types of actuators are ill-suited for the purpose 
of attitude control because of the limited range and direction of the torques they 
are capable of producing. 

Because CMGs are capable of storing and smoothly manipulating large quanti- 
ties of momentum over long periods of time, they are often favored in applications 
requiring precision pointing and the attitude control of large, long duration space- 
craft. In contrast to a reaction wheel, the small torque required along the gimbal 
axes to move the CMG rotor can produce a significantly larger torque on the vehi- 
cle. This effect is called torque amplification, and the effective gain is proportional 
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to the angular momentum stored in the CMG rotor. Also, since the CMG rotor 
spins at a constant rate, it is easier to optimize the hardware for reduced friction 
and good dynamical characteristics. In contrast to jets, CMGs do not consume fuel 
which must be periodically resupplied; rather, they use electrical power which is 
readily available in spacecraft applications. In addition, CMGs can supply a broad 
range of output torques without the quantization effects associated with jets. An- 
other undesirable characteristic of jets is that their exhaust tends to contaminate 
the environment in the vicinity of the spacecraft. 


1.1 Background 

There are two basic mechanical configurations of CMGs; the single-gimbal and the 
double-gimbal control moment gyroscope. The two types of construction are shown 
schematically in Figure 1.1. As the name suggests, the single-gimbal CMG has a 
single-gimbal drive that rotates the flywheel about an axis perpendicular to the spin 
axis. This constrains the momentum vector to rotate on a circle in R 3 . The double- 
gimbal CMG has an inner and outer gimbal drive which can be used to orient the 
momentum vector on a sphere in R 3 (in the absence of hardware limits). Single- 
gimbal CMGs tend to be lighter for a given peak torque rating and are simpler in 
construction, while control algorithms for double-gimbal CMGs tend to be easier to 
implement. 

In any particular application, control torques are generated by manipulating the 
momentum stored in a group of CMGs. The net stored momentum is given by a 
function of the form 

h = h cmg (0) (1.1) 
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Figure 1.1: Mechanical Construction of CMGs. 

where h 6 R 3 is the net CMG momentum vector in a spacecraft-fixed reference 
frame, 6 £ R m denotes the set of gimbal displacements, and h cmg : R m — > R 3 
maps the gimbal angles into the momentum space based on the types of CMGs 
involved and the geometry of the mounting configuration. This equation can be 
differentiated with respect to time in a spacecraft-fixed reference frame to yield 

r = h = LJXh + D[/i cmg (0)]0 (1.2) 

where to is the angular rate of the spacecraft, Y)[h CTng (9)} denotes the derivative of 
the mapping hcmg, and — r is the torque applied to the spacecraft due to the CMG 
gimbal motion. The torque is typically specified by an autopilot or a feedback 
control law. In general, the number of degrees-of-freedom in the CMG system, m, 
is greater than the number required to provide three-axis control of the spacecraft. 
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The problem that naturally arises is how to choose the gimbal rates, 0, to achieve a 
specified torque, r. This so-called steering problem is complicated by the inherent 
nonlinearity in these devices as well as the need to manage the redundancy in 
a CMG configuration to avoid geometric singularities. Singularities occur when 
the state of the CMG system is such that the the output torque, r, cannot be 
specified arbitrarily due to a drop in rank of the derivative in (1.2). The loss of 
rank corresponds to a loss of three-axis controllability of the spacecraft. Depending 
on the nature of the singularity, it may not be possible to rotate the CMGs out of 
the singular state without producing torques in undesirable directions [2]. 

The problem of steering single-gimbal CMGs is considerably more difficult than 
steering double-gimbal CMGs. The reason for this difficulty is that the topology 
of the momentum envelopes and singular states is more complicated for groups of 
single-gimbal CMGs [3]. For a given number of CMGs, a group of double-gimbal 
CMGs has twice the kinematic freedom which can be exploited to avoid singularities. 
Most papers on CMG steering consider only one of the two types of construction. 

To date, most approaches used in solving the CMG steering problem have in- 
volved fixing the state (u;, 9) in (1.2) and computing a set of gimbal rate commands 
that optimize some objective function subject to the constraint that a specified 
torque is produced. One exception to this method is the approach in [4] where 
gimbal displacements are computed which produce a specified angular rate change 
of the vehicle. The objective function used by the various steering algorithms is 
chosen such that its value reflects the desirability of a given configuration. At some 
point, most of these algorithms use the Moore-Penrose pseudo-inverse solution of 
(1.2) to compute a torque producing set of gimbal rates, and then add a null motion 
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[5]. A null motion is a set of gimbal rates which produces no torque, and hence 
lies in the null space of the matrix D[H cmg (9)]. The null-motion is chosen to make 
the derivative of the objective function as small - in the case of a minimization - 
as possible without exceeding physical limitations. Such an approach leads to an 
algorithm that steers the gimbals to configurations which are locally optimal for the 
given objective function. 

In [6], the objective function used by the steering algorithm is the inverse of the 
gain of a double-gimbal CMG configuration. The gain is defined to be the product 
of the singular values of 'D[H C m g {0)]i and singular configurations correspond to a 
gain of zero since they are associated with a drop in rank of this matrix. The gain 
function has been used by other authors, but has not been found to be effective due 
to the fact that it often cannot predict an impending singularity early enough to 
avoid it. It has also been demonstrated that trajectories following local optima of 
the gain can lead directly into singularities [7]. 

To attain a margin of global effectiveness, other authors have investigated steer- 
ing laws which either explicitly or implicitly locate singularities in the gimbal space 
and rotate the gimbal angles to avoid them. For a group of four single-gimbal 
CMGs, the steering law described in [7] steers to the configuration of global opti- 
mal gain for a given net momentum while avoiding tabulated singularities of the 
CMG configuration. Some disavantages of this approach are a lack of adaptabil- 
ity and the necessity of storing a very large data file. The steering laws proposed 
in [8], [5] and [4] steer the CMGs away from configurations that are known to be 
singular. The approach of [8] is to compute the sets of gimbal angles which are 
singular with respect to the requested torque vector and rotate the gimbals to avoid 
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them. In other words, gimbal angles are found for each CMG such that the CMG 
cannot instantaneously project momentum in the desired direction. The steering 
law then chooses rates that steer away from these gimbal states. The steering law 
in [5] exploits knowledge of a very specific mounting configuration to design an 
algorithm which implicitly avoids singularities. The method of [4] is probably the 
most adaptable approach in that changes in hardware and CMG configuration are 
easily accomodated, and the objective function can be changed to reflect the most 
urgent of the CMG configuration priorities. 

The algorithm proposed in [4] treats (1.2) as a constraint in a linear programming 
problem, thus avoiding some of the numerical difficulties in solving this equation 
at or near singularities. In addition, the algorithm has been used to combine the 
use of CMGs with jets to maintain vehicle control outside the torque envelope of 
the CMG system alone. Another advantage of this approach is that it explicitly 
incorporates limits on the maximum rates achievable by the gimbal drive motors. 
All other approaches that include this limitation generally impose it by scaling after 
a set of gimbal rates has been computed; thus preserving the direction chosen by 
the optimization, but violating the requirement that a given torque be produced. 
Scaling is particularly undesirable near a singular state, when the pseudo-inverse 
yields very large torque producing gimbal rates. These large rates tend to cancel 
any attempts at exploiting null-motion to avoid the singular configuration. 

In [2], a number of algorithms are investigated for the control of singie-gimbal 
CMGs. One of these methods is the singularity robust inverse method [9], which 
attempts to solve (1.2) by defining a vector consisting of the error in the torque 
produced by a choice of gimbal rates concatenated with the set of gimbal rates. A 
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solution is found by minimizing a weighted inner-product of this vector over the set 
of gimbal rates. Near singular states, this algorithm maintains reasonable gimbal 
rates by sacrificing torque in the direction associated with the singularity. 

In addition to the nonlinear behavior of the CMG system, the dynamical equa- 
tions of a rotating spacecraft are nonlinear. In particular, there is Euler coupling 
due to asymmetries in the spacecraft inertia, and coupling of the CMG momentum 
into these equations. There are also nonlinear attitude dependent torques due to 
gravity and aerodynamic effects. A number of people have linearized the equa- 
tions of motion for CMG equipped spacecraft and designed regulators for the linear 
system which optimize an integral quadratic performance measure [10], [11]. This 
approach requires assumptions that are almost never satisfied in realistic situations. 
The usual approach to such a combined spacecraft/CMG control problem is to sep- 
arate spacecraft control from CMG steering. An attitude control law generates a 
torque request that is sent to a CMG steering law, and the CMG steering law gen- 
erates a set of gimbal rate commands for the CMGs. The fact that there is a strong 
coupling between the spacecraft dynamics and the CMG motion is not exploited in 
this type of an approach. 

As noted before, spacecraft do not operate in an entirely torque-free environ- 
ment. For example, there are torques due to gravity and to aerodynamic drag in 
the very thin atmosphere, as well as forces caused by equipment or crew motion. 
In order to compensate for these disturbances, the CMGs may be required to apply 
small torques in one direction over a long period of time. Since the CMG rotor 
spins at a constant rate, each CMG in a configuration stores a fixed magnitude 
of momentum; hence, there is a maximum amount of momentum that a group of 
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CMGs can project in any given direction. As a consequence, it is necessary to ex- 
amine long-term momentum requirements of the spacecraft to insure that the CMG 
momentum capacity is not saturated. This third problem is called momentum man- 
agement. The approaches proposed for use with NASA’s space station manipulate 
the attitude dependent gravity gradient torque to control the rate at which mo- 
mentum is stored in the CMG system. Skylab used a combination of reaction jets 
and aerodynamic torques to regulate momentum build-up in the CMGs [12]. The 
solution to the momentum management problem usually employs a model to pre- 
dict the external torques acting on a spacecraft, and regulates the attitude to keep 
the resultant torque small or to produce torques which extract momentum built-up 
during previous orbits. Typically, the gravity torques can be modeled quite well, 
while aerodynamic torques are usually predicted using data from previous orbits. 
The problem of computing attitude and momentum histories to be tracked during 
an orbit has been formulated and solved as a linear quadratic optimization problem 
in [13], [14]. The linearizing assumptions may be reasonable under the assumption 
that CMGs are not included in the vehicle model, and that attitude excursions and 
rates are sufficiently small. The approach in [15] uses data from previous orbits to 
predict the momentum requirements over the next orbit due to environmental ef- 
fects and computes a sequence of attitude step commands which result in a desired 
change in momentum over the orbit. This approach has been shown to be robust 
to uncertainties in the vehicle and environmental data. 
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1.2 Thesis Objectives and Summary of Contents 


Even though CMGs have been used in a number of applications, they have not found 
widespread use. Single-gimbal CMGs have been used on balloon born platforms [7], 
in scissored pairs on early versions of NASA’s Manned Maneuvering Unit, and have 
recently been installed on the Soviet Union’s MIR space station [16]. Double-gimbal 
CMGs were used in NASA’s Skylab [17] and are proposed for use on the planned 
NASA space station [18]. The principle difficulty in using CMGs for spacecraft 
attitude control is finding reliable solutions to the steering problem. One reason for 
this difficulty is that the steering problem does not consider the long term behavior 
of the system, but uses only the current state and torque request to compute gimbal 
rates. In other words, steering laws are local in nature and do not consider the effects 
of current decisions on the future. As a result, there is no steering law which can 
guarantee that it will always keep the CMG configuration out of singular states. 

Another difficulty with using separate attitude control and CMG steering algo- 
rithms is that it is virtually impossible to predict the system behavior using locally 
computed CMG trajectories in a realistic environment. The main objective of this 
thesis is to develop an approach to the control of spacecraft with CMGs in which the 
spacecraft and the CMG actuators are treated as a combined unit, and to apply the 
approach to attitude control and momentum management of the proposed NASA 
space station. 

In this thesis, the feedback linearization theory of [19] is used in the formulation 
and solution of two important problems in spacecraft control with CMGs. The 
theory identifies a class of nonlinear dynamical systems which can be transformed 
to exhibit linear input-output behavior. The transformations do not introduce the 
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errors associated with a Taylor series approximation, and existing methodologies for 
linear systems synthesis can be applied to design closed-loop controllers, A detail 
theoretical background is presented in Appendix B. 

A number of applications using feedback linearization theory have appeared in 
the literature. In [20], [21], the theory is used to develop control laws for spacecraft 
and spacecraft with reaction wheels, respectively. Although the attitude control 
problems considered in both of these papers are quite similar, neither provides a 
realistic treatment when redundant actuators are included. In addition to the fact 
that redundancy is an important characteristic of CMG control systems, reaction 
wheel systems do not exhibit geometric singularities like CMG systems. In [22] 
and [23], the linearized dynamical systems are used in the design of autopilots for 
helicopters and aircraft, while [24] studies the control of robots. 

There is, however, a difficulty that can occur when designing control laws based 
on the linearized system dynamics. Because the transformations involve a change 
of coordinates in the state and control variables, the relationships between the lin- 
earized dynamics and the original control problem can become complicated or un- 
clear. One usually considers control problems involving a dynamical system along 
with a set of qualitative and quantitative features that characterize desirable be- 
havior. It is these features may not be easily interpretted in terms of the linearized 
system. In other words, the minimization of a quadratic performance index for 
the nonlinear system, is generally not the same as the minimization of a quadratic 
performance index with a feedback linearized dynamical constraint. In contrast to 
the developments referenced earlier, feedback linearization techniques are used in 
this thesis to simplify nonlinear dynamical models, while control and performance 
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objectives are considered in terms of the original variables. 

In applying the feedback linearization theory to the spacecraft control problem 
with CMGs, it was found that the theory is not well suited for use in some problems 
in which there are redundant actuators. The reason for the difficulty is that the 
extra degrees-of-freedom are not needed to solve the feedback linearization problem, 
and the redundancy is effectively decoupled from from the spacecraft dynamics. In 
Chapter 3, the problem of redundancy is resolved by a method that distributes 
control effort among a group of redundant actuators. This method is applied to 
the CMG steering problem to yield a class of steering algorithms in which a desired 
set of gimbal rates is computed, and then a normed approximation problem is 
solved to obtain a best match to the desired rates. The problem of designing a 
steering law is then one of determining an algorithm for computing gimbal rates 
which move the CMG rotors toward desired configurations. This is conceptually 
simpler than the conventional approach of determining a function whose derivative 
can be used to compute desired rates. Another advantage is that once a basic 
algorithm has been implemented, it is possible to add decision rules which can 
recognize and correct for known failure modes, or to design an algorithm which is 
entirely rule based. Such an expert system type of approach was employed with 
Skylab where much of the problem resolution was performed manually. In this 
thesis, a particular algorithm for determining the desired rates is proposed for the 
parallel mounting configuration to be used on the NASA space station. A number of 
examples are presented showing the behavior of the system for different formulations 
of the approximation problem. The formulation of the approximation problem 
makes use of an affine parameterization of all CMG gimbal rates that produce a 
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given acceleration. Applied to existing algorithms, this parameterization can be 
used to eliminate the torque constraint associated with the optimization, and also 
to reduce the number of decision variables by three. 

In Chapter 4, the problem of computing trajectories is addressed by formulating 
and solving an optimal momentum management problem with the CMGs included 
in the vehicle model. It is assumed that gravity gradient torques are to be balanced 
against other external torques to control the rate of momentum storage in the sys- 
tem of CMGs. All the external disturbance torques are lumped into a single term 
that is assumed known prior to solving the optimization problem or can be predicted 
using information from previous spacecraft orbits [25]. No linearizing assumptions 
are made with regard to the CMG kinematics or the dynamical equations of the ve- 
hicle. The momentum management problem is stated as the minimization of an in- 
tegral quadratic performance measure with a nonlinear dynamical constraint. This 
problem is then transformed into an unconstrained nonlinear program by feedback 
linearizing the dynamical constraint and approximating the solution of the linear 
system of equations on a reduced basis of functions. The solution of the nonlinear 
program provides control histories which produce a spacecraft attitude trajectory 
that does not require the CMGs to store excessive amounts of momentum, and 
simultaneously yields CMG trajectories that avoid singular configurations. By in- 
cluding the CMGs in the problem, it is possible to obtain a solution in which the 
CMGs avoid singular states whenever possible. Chapter 4 concludes with a con- 
ceptual description of a feedback control system architecture which will track the 
desired trajectories and reject disturbances. The design of the specific feedback 
control laws is left for future research. 
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The momentum management problem is different from the usual problem of 
rotating a spacecraft from one attitude to another. In the attitude maneuver prob- 
lem, the total momentum of the vehicle is assumed constant, or it is assumed to be 
impossible to influence the net change in momentum. The method proposed here 
for momentum management is applicable to the solution of this maneuver problem. 
The two problems are fundamentally different, however, and the application of the 
proposed method to the maneuver problem is not studied in this thesis. Since the 
NASA space station will not be required to perform large angle attitude maneu- 
vers, but will be required to hold a specified attitude for long periods of time, the 
momentum management problem is the appropriate trajectory planning problem 
to study. Results presented in the following chapters are derived for general space- 
craft/CMG configurations, but the specific examples all use vehicle data and CMG 
configurations for the planned NASA space station. 

Some results in this thesis are applicable to problems other than the control of 
spacecraft with CMGs. In particular, there is a strong analogy between the steering 
problem for CMGs and the inverse kinematics problem encountered in robotics. The 
connection is made in [2] for single-gimbal CMGs and planar, serial-link robots, 
the relationship between singularities in the two systems is made quite nicely in 
this thesis. The robotic analogy for a double-gimbal CMG is a robot joint with 
two degrees-of-freedom. Most steering laws can be interpreted as resolved motion 
methods when applied to inverse kinematic problems involving robots [26]. The 
attitude maneuver problem discussed in the preceding paragraph corresponds to the 
problem in robotics of computing optimal inverse kinematic trajectories which take 
a robot end-effector from an initial orientation to a specified final orientation. The 
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method of Chapter 4 is applicable to this problem. In fact, a recent paper describes 
an approach similar to that taken in Chapter 4 for the robotic problem [27]. In this 
paper the authors use penalty function methods to satisfy the constraints imposed 
by the dynamical equations on the nonlinear program. The approach of Chapter 4 
could be used to satisfy the constraint exactly and to state the problem in terms of 
a significantly reduced number of variables. 

In general, the approach to optimization outlined in Chapter 4 can be applied 
to all nonlinear systems which are feedback linearizable as defined in Appendix B. 
Broad classes of models for aircraft and robotic systems are known to satisfy these 
conditions. For systems which cannot be linearized, provisional methods are out- 
lined for the solution of the resulting constrained nonlinear programs. 

Before developing the feedback linearization theory and applying it to the two 
general problems considered in the thesis, equations of motion are developed in 
Chapter 2. Chapter 2 also introduces the methods of describing rotation with 
quaternions and presents necessary elements of quaternion algebra. Equations of 
motion are developed which describe the attitude dynamics of a spacecraft free 
to rotate without the influence of feedback torques. Equations are subsequently 
developed in which the gravity gradient torque is included. The chapter includes 
a brief discussion of the feedback linearization problem and a method of solving 
a simplified version of the problem. The chapter concludes with a derivation of 
transformations that linearize the two problems. The first torque-free system of 
equations is used in Chapter 3 to study the steering problem. The augmented set 
is applied in Chapter 4 in which an optimal momentum management problem is 
stated, and a technique for obtaining solutions is developed. 
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Chapter 2 

Derivation of Equations of Motion 


In this chapter, equations of motion are derived for a rigid, rotating spacecraft 
equipped with CMGs. The attitude equations are written in terms of the quaternion 
relating a reference frame, fixed to the spacecraft, to an inertial frame. The next 
section contains an overview of relevant facts regarding quaternions and quaternion 
algebra. Kinematic models are then derived for single- and double-gimbal CMGs, 
and appropriate reference frames are defined for expressing the momentum of a 
CMG in a spacecraft fixed coordinate system. The reference frames depend on the 
adopted mounting configuration and the type of CMGs employed. In the third 
section, equations are derived that describe the attitude dynamics of a spacecraft 
undergoing torque-free rotation and subject to an internal transfer of momentum. 
The section concludes with the development of a system of first-order differential 
equations that fully describe the dynamics of the combined spacecraft /CMG system. 
In the next section, equations of motion are derived that include the effects of gravity 
gradient torques acting with a general external torque applied to the vehicle. In 
the final section, a brief introduction to the linearization problem is given, and 
linearizing transformations are developed for the two spacecraft models that are 
derived in this chapter. 
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2.1 Quaternions and Quaternion Algebra 


There are a number of ways to characterize finite rotations of rigid bodies - or ref- 
erence frames - in orthogonal coordinate systems. For example, Euler angles and 
direction cosines are common choices of coordinates in many applications. It is 
useful to have a representation for rotation that has an algebraic operation corre- 
sponding to successive rotations of a reference frame. Quaternions were invented 
by Sir William Rowan Hamilton in 1843 to extend ordinary 3-dimensional vector 
algebra to include operations of multiplication and division. One important use of 
quaternions is to transform vectors expressed in one reference frame into the ba- 
sis of another reference frame. The algebra-like structure of quaternion arithmetic 
makes it easy to perform successive transformations between reference frames. A 
brief introduction to the relevant aspects of quaternion algebra is presented in this 
section; see [28] for additional details. 


2.1.1 Definition and Elementary Operations 

A theorem proven by Euler states that any sequence of rotations of a rigid body 
- or of a coordinate frame - is equivalent to a single rotation about some axis. 
Quaternions are essentially mathematical tools that implement this theorem. Let 
A and B denote a pair of right-handed reference frames, and assume that B has 
been rotated with respect to A through an angle 6 about an axis defined by the 
unit vector u expressed in A . The quaternion of reference frame B with respect to 
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A is defined to be the ordered quadruple 

cos(0/2) 
usin(0/2) 

Quaternions are usually regarded as a quadruple composed of a scalar component, 
q 0 , and a vector component, q, consisting of the last three elements of the quaternion. 
Definition 2.1 is more restrictive than typically employed in general discussions 
[28]. This particular definition is useful when quaternions are only being used to 
transform vectors between reference frames. 

Quaternion addition and subtraction are the ordinary operations for ordered 
n-tuples of numbers. Quaternion multiplication is defined by 

Po qo-p-q 

P o9 + 9oP + P x 9 

where • and x denote the ordinary vector dot and cross products, respectively. The 
multiplication operation defined by (2.2) has the following interpretation. If q% 
denotes the quaternion of frame B with respect to A, and q% denotes the quaternion 
of C with respect to B, then their product 

q c A = p c b q B A (2-3) 
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is the quaternion that defines the orientation of frame C with respect to frame 
A. The definition of quaternion multiplication can be expressed in vector-matrix 
notation to yield 


■ 






P 0 

“Pi 

“P 2 

~P 3 


9o 

Pi 

P 0 

~P 3 

P 2 

< 

9i 

> 

P 2 

P 3 

Po 

“Pi 


92 

P3 

“P 2 

Pi 

P° 


. 93 , 


(2.4) 


as the definition of quaternion multiplication. The conjugate of a quaternion q is 
defined to be 


9 


* 


A 



(2.5) 


It can be shown that quaternions which perform transformations between orthogo- 
nal reference frames have the property that 


9a 



( 2 . 6 ) 


Using the definition of quaternion multiplication, it can be shown that the quater- 
nion i is the multiplicative identity for quaternion algebra. The implication of (2.6) 
is that the conjugate of a quaternion defined by (2.1) is in fact its inverse. This is 
analogous to a similar result for transformation matrices. 
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2.1.2 Coordinate Transformations with Quaternions 


Quaternions can be used to transform vectors expressed in one reference frame into 
another coordinate frame. A vector x A which is coordinated in a reference frame 
A can be tranformed into another reference frame B by the sequence of quaternion 
products 


x 


B 


„B* A „B 

9a x 9 a 


(2.7) 


where the vectors x A and x B are expressed as quaternions whose scalar parts are 
zero. It is the simplicity of this expression which motivates the use of (2.1) as the 
working definition of a quaternion. The effect of this transformation is to conically 
rotate the vector x A about an axis defined by g, or equivalently u, through an angle 
0 . Another way of viewing this transformation, is to consider what an observer 
in reference frame B would see as his reference frame was rotated from a position 
initially aligned with A . To this observer, it would appear that the vector x A was 
being rotated about the axis — u. From a computational point of view, the most 
efficient implementation of the coordinate transformation (2.7) of a vector using 
quaternions is 


X B = (2?o 2 - 1) X A + 2 (q • X A ) q - 2q 0 (qx * x ) 


( 2 . 8 ) 


Successive transformations of a vector between several reference frames can be 
computed by repeated application of (2.7) 


c C* ( B* a B \ c 
x = 9b {9a x 9a) 9b 


(2.9) 
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Since it can be shown that ( pq )* = q*p *, this equation can be rewritten as 


& * A (dti) 


„C* A C 

9a x 9a 


This last equation reiterates the connection between multiplication of quaternions 
and rotations of orthogonal reference frames. The relationship between quaternions 
and the direction cosine matrix associated with a rigid rotation is given by 

2 (?o + 9 i) — 1 2(5152 + 9093) 2(5153 — 9092) 

C{q) = 2(qiq 2 -qo93) 2(5 q + 5!) — 1 2(5293 + 9o9i) C 2 - 11 ) 

2(5153 + 5052) 2(5253 - 5o9i) 2(5^ + 53) - 1 

where the fact that 5o + 5i + 52 + 93 = 1 has been used to simplify the diagonal ele- 
ments. Equivalently, a quaternion can be extracted from a direction cosine matrix, 
C = [cy]. The diagonal elements of (2.11) yield the system of equations 

50 = 1 + Trace(C) 

5? = 1 + 2cn — Trace(C) 

( 2 . 12 ) 

q\ = 1 + 2 c 2 2 — Trace(C) 

51 = 1 + 2 c 3 3 — Trace(C) 

and the off-diagonal elements yield 


PoPl 

— c 23 

— C32 

P2P3 

= C 2 3 + C32 

P0P2 

— C 31 

— c 13 

P3P 1 

= c 31 + C 13 

P0P3 

“ c 12 

— C 21 

PlP2 

= c 12 + C 21 


( 2 . 13 ) 
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Choosing a nonzero q{ from (2.12) and taking its positive square root makes is 
possible to solve three of the equations in (2.13) for the remaining elements of the 
quaternion. 


2.1.3 Derivative of a Quaternion 


In the following sections, it will be necessary to compute the time rate of change of a 
quaternion due to the relative rotation of two reference frames. For example, if the 
quaternion transforms vectors in inertial coordinates into a spacecraft fixed frame, 
the quaternion will change with time as the spacecraft rotates. In this section, the 
derivative of a quaternion, q B , is computed assuming reference frame B is rotating 
at a relative angular rate w B which has been coordinated in the reference frame 
B. 


In order to determine the derivative of a quaternion, it is necessary to evaluate 
the limit 

_ a /\n ^ ( i 4 - ( i\ 

(2.14) 


&(t) = lim «4(< + A«)-^( j) 

' At— ^0 A t At-0 At 


After an interval At , the reference frame B will have rotated through an angle 
u; 5 |A< to a new orientation B f about a unit vector Q B aligned with u > B . The 
quaternion which performs this transformation is 



cos(|a> B | At/2 ) 

► 


Q b sin(|o; B At/2) 


(2.15) 
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and the quaternion q B (t + A t) is given by 


q%(i + At) = q% q% 

The limit (2.14) can now be evaluated from 


(2.16) 


= lim 


d'd - d 


A t 

= la h 


= lim — < 
Ae->o At 


cos(|u; B At/2) — 1 




Q b sin( u> B At/2) 


1 

= 2< 


a ; 


B 


9a 


(2.17) 


where i denotes the unit quaternion, (2.6). 

It is also of interest to determine the relationship between the relative accel- 
erations of two reference frames and the second time derivative of the quaternion 
relating them. Expanding (2.17) using the definition of quaternion multiplication, 
(2.2), and differentiating yields 


9 n Ji 


l_d 

2 dt * 


-u ) B • q 


ui B q 0 + u) B X q 


1 

2 


* B - 

1 

+ - 

B — 1 

— u> * q 

i 

-q | 

/ 

Cj b qo + x q 

i. j 

2 

U) B q'o +w B X q J 
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I 


0 


0 


9 


LO 


B 




LO 



(2.18) 


where the subscript and superscript on the quaternion have been dropped to simplify 
the notation. 


2.2 CMG Kinematics 

In order to derive dynamical equations for a spacecraft equipped with CMGs, it 
is necessary to develop a kinematical model for CMGs. In this section, the equa- 
tions relating CMG gimbal angles and CMG momentum are derived. Since the 
mechanical construction of single- and double-gimbal CMGs are quite different, the 
relationships are derived separately. In each case, a reference configuration is de- 
fined which is used as a coordinate system for the CMG rotor [29]. The kinematical 
constraints are then used to express the CMG momentum vector in this reference 
frame as a function of the gimbal angles. 

Figure 2.1a shows a schematic drawing of a single-gimbal CMG. The unit vector 
a lies along the axis about which the CMG is rotated by the gimbal drive, and 
this vector is fixed with respect to the spacecraft. This axis is called the gimbal 
axis and the displacement of the rotor about this axis is denoted by 9 . The unit 
vector h° corresponds to the orientation of of the CMG momentum vector at the 
reference configuration, 0 = 0. Together with the vector a x /i°, these two vectors 
form a right orthogonal reference frame. The orientation of the CMG rotor in this 
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reference frame is 

h(0) = h° cos 9 + {S’ x h°) sin 9 (2.19) 

and the momentum of the CMG in the basis (?, /i°, a x h° ) is given by 

h sg (9) = h mag ^hP cos 9 + (<? x h°)sin0] (2.20) 

where h mag is the magnitude of the momentum stored in the rotor. 

The equations describing the kinematics of a double-gimbal CMG are derived 
using the same principles as above. The difference is that the rotor of a double- 
gimbal CMG is rotated about two axes and the orientation of the momentum vector 
depends on two parameters; the outer gimbal displacement, 9 0 and the inner gimbal 
displacement, 0;. Figure 2.2a shows a schematic of an Euler mounted double-gimbal 
CMG. In this figure a 0 is a unit vector along the outer gimbal axis. This axis is 
fixed with respect to the spacecraft. The unit vector lies along the reference inner 
gimbal axis corresponding to 9 0 = 0, and the third vector in the triad is defined 
to be orthogonal to a Q and a 0 (i.e. h° = h(9 0 = 0,0; = 0) = a 0 X 5^ ). The 
CMG momentum vector can now be expressed in the basis (<r 0 , <?*, /i°) shown in 
Figure 2.2b. 

Holding the inner gimbal fixed at zero and rotating the CMG about the outer 
gimbal axis yields the expression 

h(0 o ,O) = h°cos0 o — a{ sin0 o (2.21) 

for the the unit vector along the rotor axis. Rotating the inner gimbal about the 
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a.) Schematic. 



h{o oy o y h{0 o ,di) 


b.) Reference frame. 


Figure 2.2: Schematic Drawing and Reference Frame for a Double-gimbal CMG 
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rotated inner gimbal axis, yields the expression 


h(9 0 , 0,) = h{9„ , 0) cos 0,- + <r 0 sin 0,- (2.22) 

Substituting (2.21) into (2.22) and scaling by the magnitude of the rotor momentum 
yields 

hdg{9o, 0.) = hmag [(h° cos 9 a - sin 0 O ) cos 0,- + <r 0 sin 0,] (2.23) 

as the final expression for the momentum of a double-gimbal CMG. 

Equations (2.20) and (2.23) express the momentum vector of a single-gimbal 
and a double-gimbal CMG, respectively. The momentum vectors are expressed in 
reference frames which are defined by the mounting of the individual CMGs. In 
general, these reference frames will be coordinatized in a common frame which is 
fixed with respect to a spaceraft. Given a group of CMGs, the equations derived 
in this section can be used to express the net momentum of the CMG system in a 
spacecraft-fixed reference frame as a sum over the set of CMGs: 

N cmg 

W») = E M«<) (2.24) 

1=1 

where h{ £ R 3 are the momentum vectors of the individual CMGs, 6 £ R m is the 
set of gimbal angles, and d{ denotes the gimbal angle(s) associated with CMG i. 
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2.3 Equations of Motion for a Rigid Body: No External 
Torques 

In this section, the equations of motion are derived for an isolated rigid spacecraft 
rotating in an inertial reference frame. The assumption of isolation implies that 
there are no external torques acting on the spacecraft and that the total momentum 
of the spacecraft is constant when expressed in inertial coordinates. Anticipating 
the need to include CMGs in the dynamical equations, the total momentum of the 
spacecraft is regarded as the sum of the momentum associated with the angular 
rate of the spacecraft and the momentum stored in a system of CMGs on board the 
vehicle. The equations are derived by a straight-forward application of conservation 
of momentum. 

At this point it is useful to make a qualification on the definition of quaternion 
multiplication. Throughout this thesis, it will be necessary to compute quaternion 
products of only two types of quaternions: Quaternions associated with transfor- 
mations between orthogonal reference frames, and quaternions associated with 3- 
vectors to be transformed. The scalar part of a quaternion representing a 3-vector 
is always zero. Since there cannot be more than three independent quantities which 
characterize rotations in a three-dimensional space, the four elements of the quater- 
nion cannot be independent. A consequence of definition (2.1) is that quaternions 
which perform transformations between orthogonal reference frames have a Eu- 
clidean norm of one; given any three elements of such a quaternion, it is possible to 
compute the fourth. In particular, the expression 

?o = . 1 - £ q} (2.25) 

\ *=1 
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can be used to compute the scalar part of a quaternion from its vector part. From 
now on, it will be assumed that all quaternions are represented by triples of num- 
bers; hence, the definition of quaternion multiplication (2.2) must be modified to 
use either zero or (2.25) when the scalar part of the quaternion is needed. With or 
without subscripts, the symbol q is reserved exclusively for quaternions that per- 
form orthogonal transformations, and quaternions associated with 3-vectors will be 
0 

denoted by >. The notation x will also be used when the interpretation is clear 
x 

s * 

from the context. With these definitions, it is possible to derive a minimal system 
of first order equations to describe the system. 

The attitude of the spacecraft is characterized by a quaternion which defines the 
relative orientation of a spacecraft-fixed reference frame, denoted J9, with respect 
to an inertial frame, gj . The total angular momentum of the spacecraft in inertial 
coordinates, /i 1 , can be transformed into the basis of the spacecraft-fixed frame 
using (2.7) to yield 

h B = q* r h 1 qj (2.26) 

Note that while h 1 is a constant because there are no external torques, h B changes 
in time as the spacecraft rotates due to the evolution of the quaternion, g/. The 
momentum, h B , can be expressed in the spacecraft reference frame as 

h B = h CTng (6) + h 8pacecra f t = h CTng (0 ) + Iuj b (2.27) 

where h cmg (9 ) denotes the net momentum stored in the system of CMGs, co denotes 
the spacecraft’s angular rate, and I denotes the inertia tensor of the vehicle. Both 
h CTng (0)^ and u; are expressed in the body fixed coordinate frame B . The angular rate 
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of the spacecraft can be determined from the total momentum and the momentum 
stored in the CMGs by solving (2.27) for u> 

u, B = I" 1 (h B - h^O)) = I- 1 (q! h 1 qj - h-tf)) (2.28) 

This expression can be substituted into (2.17) to obtain an expression for the deriva- 
tive of the attitude quaternion in terms of the constant inertia matrix, the constant 
inertial angular momentum, and the attitude quaternion of the spacecraft. 



(2.29) 


In order to complete the derivation of the equations for rotation in the absence 
of external torques, it is necessary to consider the control or input variable. CMG 
steering algorithms usually choose the set of gimbal rates as the control variable. 
Since the instantaneous torque produced by a group of CMGs is a function of the 
gimbal rates only, the gimbal rates are the only choice available when the steering 
law is required to produce specified torques. The gimbal rates are also useful when 
it is desirable to include physical limits on the individual gimbal drives. ; In some 
circumstances, it may be more realistic to consider torque limits on the gimbal drive 
motors, especially when double-gimbal CMGs are being used. With single-gimbal 
CMGs, the torque applied to the vehicle is orthogonal to the axis of rotation of the 
rotor. With double-gimbal CMGs, motion of one gimbal can produce a torque along 
the axis of the other. The need to back-drive the gimbals can impose restrictions 
that can not be accounted for without developing a higher-order model. In this 
discussion, it is assumed that that the gimbal rates - as expressed in the individual 
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CMG frames - can be commanded as the control variables. This assumption is 
consistent with conventional assumptions and does not result from a limitation 
imposed by the feedback linearization theory that will be applied to the problem. 
It may be desirable in future efforts to examine the effects resulting from including 
the dynamics of the gimbal drive motors in the model. These phenomena are not 
significant at the low torque and bandwidth levels encountered in this study and 
are not incorporated. 

Under the stated assumptions, the dynamics of the combined spacecraft /CMG 
system are described by 

9 = u 



(2.30) 


where 9,u 6 R m . These ra + 3 equations describe the evolution of the state variables 
for this system of equations: The three elements of the attitude quaternion, q/, and 
the m gimbal angles, 6. 


2.4 Equations of Motion for a Rigid Body: External Torques 

In this section, we drop the assumption that the spacecraft is an isolated rigid 
body. The external torques that act on a space vehicle are generally small; however, 
the effects of these torques can be significant over long periods of time. The two 
principal sources of torques on orbiting spacecraft are due to gravity-gradient and 
aerodynamic effects. The attitude of the spacecraft will be specified with respect 
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Figure 2.3: Orientation of LVLH Reference Frame 

to a local vertical /local horizontal (LVLH) reference frame in a circular orbit about 
the Earth. The origin of this reference frame is fixed to the center- of- mass of the 
spacecraft. The LVLH reference frame is defined as shown in Figure 2.3. The 
x-component of the triad is aligned along the velocity vector, the y-component is 
aligned with the angular velocity vector of the orbit, u; 0 , and the ^-component is 
aligned with the radius vector of the orbit. With this definition, the LVLH reference 
frame defines a plane which is locally horizontal with respect to the Earth and a 
direction which is locally vertical The orientation of the body frame with respect 
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to the LVLH reference frame will be denoted by the quaternion q 

In order to compute the derivative of the attitude quaternion, it is necessary to 
determine the angular velocity of the vehicle relative to the LVLH frame. Assuming 
that the net angular momentum is known, (2.28) can be used to determine the 
inertial angular rate of the vehicle in body coordinates, u* B . The angular rate of 
the spacecraft with respect to the LVLH frame, a; L , can be computed in LVLH 
coordinates by expressing d b in the LVLH basis and subtracting the angular rate 
of the orbit, u) Q : 


UJ L = q L U} B q* L - lo 0 

= q L {I" 1 (ql h L q L - h cmg ) } q* L - (2.31) 

where h L denotes the net momentum expressed in the LVLH basis. 

The gravity-gradient torque arises because the forces on the orbiting vehicle 
due to gravity and centripetal acceleration vary as different functions of the orbital 
radius, R 0 . In particular, the centripetal force is proportional to while the 
gravity force is proportional to R ~ 2 . If a spacecraft is orbiting the Earth in a circular 
orbit, the radial force due to centripetal acceleration and gravity exactly balance at 
the center-of-mass. For a small element of mass orbiting at a distance greater then 
the center-of-mass, the centripetal force is larger than the gravity force and there 
is a net force in the outward direction. The opposite is true for an element of mass 
orbiting at a lower altitude. For a long cylindrical object in a circular orbit, the 
gravity gradient torque will be zero if the longitudinal axis is aligned with the local 
vertical. The expression for the gravity gradient torque in a body-fixed reference 
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frame in orbit about a spherical primary body (i.e. an ideal Earth) is [30] 

T 9B = ( 2 - 32 ) 

where /x = 3.986012 x 10 14 m/ sec. is a constant and 6 denotes a unit vector pointing 
to the center-of-mass of the primary. 

Since external torques act on the spacecraft, it is no longer valid to assume that 
angular momentum is a constant. Newton’s Third Law can be applied to obtain an 
expression for the rate-of-change of momentum caused by the external torques. The 
gravity gradient torque (2.32) can be expressed in LVLH coordinates by substituting 


o h 


0 

< 0 
-1 


and I L = C T IC 


(2.33) 


where d L denotes the vector to the center of the orbit, I L denotes the vehicle inertia 
tensor in the LVLH basis, and C is the transformation matrix (2.11). Using this 
substitution, the net external torque acting on the vehicle in the LVLH basis is 


r 


L 


3 fi 

K 


-i L 

1 23 

I L 

J 13 

0 


“ 1 “ QL T"ext 


(2.34) 


here r ext denotes all other external torques on the vehicle (e.g. aerodynamic, crew 
push-off, etc.) and is expressed in body coordinates. 

As defined in the previous section, the gimbal rates will be considered the input 
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variable. A system of first order equations can now be derived using (2.31) with 
(2.17), and using (2.34) with the fact that the apparent rate of the inertial reference 
frame with respect to the LVLH frame is —w 0 , to yield: 

0 = u 



0 

I Q L 

!~ l ['ll <1L - hcmg{0)] j 

These ten equations describe the time-evolution of the attitude quaternion, the 
net momentum in LVLH, /i L , and the gimbal angles, 9. 
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2.5 Linearization of the Two Spacecraft Models 


Beginning with the work of Krener [31], many researchers have investigated the 
conditions under which nonlinear dynamical systems can be transformed locally 
into linear, controllable systems. The transformations are local in the sense that 
they exist at all states sufficiently close to a nominal state. Depending on the 
specific types of transformations being allowed, different classes of systems can be 
linearized. Unlike Taylor series approximations, the linearized model is an exact 
representation of the nonlinear system. 

Krener introduced and solved the linearization problem when the allowed trans- 
formation was only a local change of coordinates in the state space [31]. This 
state-space linearization problem was generalized to allow a linear change of input 
coordinates in [32]. The authors of [19] found a class of systems which can be 
transformed to a linear system using state dependent feedback, linear coordinate 
changes in the input, and nonlinear state coordinate changes. Figure 2.4 illustrates 
how the feedback and coordinate changes define a transformation which linearizes 
the nonlinear dynamics of a system. The figure shows that the transformation is 
applied externally to a nonlinear system; in theory, these transformations could 
be implemented directly in hardware and/or software for application to a physical 
system. 

Since the foundations of this problem were established, a number of authors have 
extended it to include a wide variety of generalizations. The problem of transforming 
nonlinear systems with specified outputs using state coordinate changes is solved in 
[33], and a class of systems which are partially linearizable is studied in [34]. The 
problem of globally linearizing scalar input systems is solved in [35]. These ideas 
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u = a (x) + p(x) v 

x = f(x) + G(x) u -< ► z = A z + B v 

z = F(x) 


Figure 2.4: Implementation of a Linearizing Feedback Transformation. 
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have also been extended to the dual problem of transforming nonlinear systems 
with output measurement into systems for which linear, spectrally assignable state 
estimators can be built [36]. 

Rather than treat the general feedback linearization problem here, a simpler ver- 
sion of the problem will be defined and solved. The general problem along with a 
detailed background discussion can be found in Appendix B. The problem discussed 
here is the feedback linearization problem with output specified. In this formulation 
of the linearization problem, all of the state variables of the linearized system are 
effectively specified at the outset. If the system is linearizable in terms of these 
coordinates, then the constructive procedure for obtaining the linearizing transfor- 
mation is quite simple. The question of existence of such a set of state variables 
is far more delicate, and is considered in detail in Appendix B. The procedure for 
solving the problem discussed here is then applied to linearize the two dynamical 
spacecraft models that have been derived in this chapter. 

2.5.1 Feedback Linearization with Output Specified 

In this section, linearizing transformations are obtained for nonlinear systems of the 
form 


x = f(x) + Yj 9i{ x ) u i = f{ x ) + G(x) u (2.36) 

t=l 

y = h{x) (2.37) 


where the nonlinear, vector valued output equation y = h(x) has been added. It 
is desired to transform this equation so that it exhibits the linear behavior of the 
system shown in Figure 2.5. The system of decoupled strings of integrators shown 
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m-1 




(m-l)2 (m-l)l 



Figure 2.5: Brunovsy Form for the Linearization with Output. 
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in the figure is known as a Brunovsky cononical form . Suffice it to say that all con- 
trollable linear systems can be brought to this form, and if the linearizaton problem 
is solvable at all, the nonlinear system can be made to behave like the system in 
Figure 2.5 (refer to the appendix for details). The indices K; are the Kronecker 
indices of the linear system [37], and characterize the number of integrators in each 
row of the canonical form. 

It is assumed that the output functions = h{{x ) in (2.37) are to correspond 
the outputs of the strings of decoupled integrators in the Brunovsky form. If there 
is a transformation that maps the dynamics of the nonlinear system to the system of 
Figure 2.5, then the functions that map x to the state variables of the linear system 
are in a sense time derivatives of the functions /i,(x). Assuming that the linearization 
problem is solvable with these outputs, then the coordinate transformation and 
the feedback that linearize the system can be obtained by differentiation. In the 
following paragraphs, the coordinate tranformation z = F(x ), and the feedback 
u = a(x) + (3{x)v that transform the nonlinear equations (2.36) and (2.37) to the 
canonical form shown in Figure 2.5 are derived. 

In order to simplify some of the following equations, the notation L w <f>(x) is used 
to denote the Lie derivative of the scalar valued function </> : R n — > R with respect 
to the vector valued function tu(x). The Lie derivative is computed according to 
the rule 

L w <f>(x) = Df>(x)w(x) = ^2 Wi(x ) (2.38) 

»=i 0/1 

and is simply the directional derivative of the function (j)(x ), in the direction tn(x), 
at the point x . 

Suppose the component function hi(x) in (2.37) corresponds to the state variable 
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t/ix in Figure 2.5. One component of the coordinate transformation F(x) is clearly 
the function hx(®). Since y 12 is the time derivative of y u , the function that maps 
to yi 2 in Figure 2.5 must be given by 

Vi 2 = Vu - Dhx(®) X = Li hi(x) (2.39) 

Substituting the expression for x given by (2.36) yields 

Lji/ii(®) = Lj+cuh^x) = Lfhi(x) + Lc tl /i 1 (®) (2.40) 

It should be clear from the figure that yu is functionally independent of the input 
to the linear system, v. Since u and v are related algebraically, it must also be true 
that t/ 12 is independent of u, implying Lg u h\(x) = 0 for all choices of u. This leads 
to the conclusion that y i2 = Lfhi(x), By succesive application of the Lie derivative 
operator L f to the function hx(x), the value of each of the state variables in the 
first set of integrators can be obtained according to the rules 

y u (x) = hx(x) 

Vu{x) = yn = L f hi(x) 

Vi3{x) = y \2 =LfL f hi(x) ( 2 - 41 ) 

2/l«i ( ^ ) = ili(Ki-l) —h j hl(^) 

In a similar fashion, the functions that map the variable x to each of the yij in 
Figure 2.5 can be obtained. By ordering these functions in a vector of the form 
z — • • ■ »2/m,/t m ] j the mapping z = F(x) can be constructed. 
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The state transformation was constructed by evaluating Lie derivatives with 
respect to x of the functions hi(x ) up to order — 1. Inspection of Figure 2.5 
indicates that evaluating the Lie derivative of the left-most state variable, y t(ti of 
the t th string of integrators will yield an expression for the i th input to the linear 
system, Vi = The manipulations yield a set of functions which are not 

independent of u 

V{ = (®) = Tjf+Q u Tjf hj(®) 

= L fV?- x hi{x) + L Gu LJ- l hi{x) 

i (2-42) 

= i) i h i {x) + 'L% l i JgiUi iy- l hi(x) 

m r m /}?/• 

= l ;%(*) + £ 

i= i U=i OXk 

where denotes the k th element of the column vector gi(x). The m expressions 
for the individual components can be combined to obtain a single vector-matrix 
equation for the control vector v 

v = w(x) + A(x) u 

where 

l 7M*) | 

i and j4(b) = [aij(x)] = £ IT 51 £*>(*) u (2.44) 

U=l OXk 

L/ m ymi(x) _ 

If the linearization problem is indeed solvable, then the matrix A(x) is nonsingular 



(2.43) 
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and (2.43) can be solved for the u which produces a desired v. Defining 


u 


a(x) + (3(x)v where 


a{x) 


0{x) 


-A-'(x) 


Vf'V u(x) 


A-'(x) 


(2.45) 


the feedback linearization problem is solved by the state coordinate change F(x), 
the state feedback a(x), and the linear change of coordinate in the input space 0(x). 

From the preceeding discussion, it is clear that a solution to the feedback lin- 
earization problem can be obtained if the structure of the Brunovsky form in Fig- 
ure 2.5 is known (characterized by its Kronecker indices), and a set of m functions 
y,i are known which are the outputs of the strings of integrators in the figure. In the 
following two subsections, linearizing transformations are derived for the two mod- 
els developed in the previous chapter. In both subsections, the Kronecker indices 
are given along with a set of functions which solve the linearization problem. 
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2.5.2 Linearization of the Spacecraft Model without External Torques 


The spacecraft model (2.30) can be put in the standard form for the linearization 
problem by letting 


/(*) = 



0 



c \ 
Ou 

1 

- 

2 

0 

f > 

0 

/-■ [qjh'q, - />„,,(«)] J 

1 qi 

) 

and 5 j(x) = 

&mi 

0 

0 

V 0 / 


(2.46) 


with the nonlinear system to be linearized described by the vector fields f,gi,...,g m . 
The first m elements of the state vector x are the gimbal angles, and the last 
three elements of x are the vector components of the attitude quaternion. Direct 
computation shows that the Kronecker indices of the spacecraft model (2.30) are 


Ki 


2 for i = 1,2,3 

< 

1 for i = 4, . . . ,m 


(2.47) 


The Brunovsky form is a set of 3 double-integrators and m — 3 single-integrators. 
The double-integrators correspond to the rotational dynamics of the spacecraft, 
and the single integrators are related to the rotation of gimbals. While the three 
double-integrators should have been expected from the physics, the number of single 
integrators may seem unusual. The reason for this occurence was alluded to in 
the introductory paragraphs. Once an attitude variable is chosen, it is necessary 
to include its time-derivative as a state variable. The derivative is algebraically 
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related to the angular rate of the vehicle, and therefore, to its momentum. The 
fact that the total inertial momentum is fixed necessarily fixes the net momentum 
of the system of CMGs. Constraining the momentum of the CMG system implies 
a loss of 3 degrees-of-freedom; hence, the number of gimbal-like control variables in 
the linearized system is m — 3. 

One choice of functions y»i (a:) that solves the linearization problem is 

yn(s) = 9i 
x) = q 2 

2/31 (*) = ?3 

2/«i(z) = Oi for i e I* 

where I* denotes an index set consisting of m — 3 integers from 1 through m. The 
three indices which are dropped must correspond to three gimbals which together 
provide three axis torque control of the spacecraft. The remaining state variables 
are computed by Lie differentiation 

Vi2{x) = 


for £ = 1,2, 3. The linear change of coordinates in the input space is computed from 


y.i(z) 

I'/yii(®) 

E fj( x ) 


j=m + 1 


dyn 

dxj 


(2.49) 



(2.48) 
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(2.43) by defining 


fhtM 


ddi 

00m 

d V22 

02/22 

00! 

00m 

§ 2/3 2 

01/32 

00! 

00m 

••• 

8i\m 


(2.50) 


L ••• 

where the notation ij denotes the j lh element of the index set I*. A second set of 
Lie derivatives is computed according to 


Wi(x) = L)yn{x) 




n 

= E 

j=m+l 


r 9hjyi2 
Sj dxj 


(2.51) 


for i = 1,2,3. The linearizing feedback and input coordinate change are given by 


a(z) = —A l (x)w(x) and = A l (x) (2.52) 

The state variables of the linearized system are the vector component of the 
attitude quaternion, the derivative of the attitude quaternion, and a subset of the 
available gimbal angles. The control variables associated with the linear system are 
the second derivative of the attitude quaternion and a subset of the available gimbal 
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rates. In terms of the notation of the previous section, ( z,v ) are given by 


A 
Z = 


( \ 

q 




and v = 


'i' 

\°‘ / 


(2.53) 


This choice of coordinates is by no means unique; any independent set of attitude 
variables would work. In most attitude control problems, it is customary to think 
of applying torques to a vehicle to cause rotation. The fact that r = ICo B to- 
gether with (2.18) provide the algebraic connection between commanding torques 
or accelerations, and commanding a quaternion second derivative for the purpose 
of control. 


2.5.3 Linearization of the Spacecraft Model with External Torques 

The spacecraft model that includes gravity gradient torques can be put into the 
standard form by letting 


/ 


/(*) = 


1 

2< 


0 


\ 


. . l L fy* JL 

3 03^23 

Kq 

3 /* jl 

D3 i 13 

K 0 
0 

[ I~ l [q L hL QL - h C mg(0)] J 


} qi 


and gi(x ) = 


l , \ 

0 li 


Smi 


0 

0 


\ 0 / 


(2.54) 
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with the nonlinear system to be linearized described by the vector fields /, <7i, . . . , g m - 
The first m elements of the state vector x are the gimbal angles, the next three 
elements are the components of the total momentum in LVLH, and the last three 
elements of x are the vector components of the attitude quaternion. The term 
associated with the external torque, qL^ext^Li has not been included in (2.54) because 
the linearization problem is solved for a set of coordinates in which this term appears 
linearly. In this way, the external torque appears in the linearized system as a 
disturbance, rather than being cancelled automatically by the feedback associated 
with the linearization. 

The Kronecker indices for the spacecraft described by (2.54) are 


4 


for i = l 




3 for i = 2 
2 for i = 3 


1 for i = 4, . . . , m 


(2.55) 


The Kronecker indices of this system clearly indicate that the relationships between 
the state variables are more complicated for the case in which gravity gradient 
torques are included. As in the previous case, three gimbals must be used to provide 
attitude control of the spacecraft; therefore, only m — 3 of the available gimbal angles 
can be associated with state variables in the linearized system. The triple-integrator 
and the first three integrators associated with k,\ appear because the rate of change 
of momentum is related to the attitude. In other words, the momentum is given 
by an integral of a function of the attitude. There is a double-integrator in the 
linearized system because no attitude-dependent gravity gradient torque acts about 
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the z direction in LVLH axis. 

Before proceding with the linearization, it is usefull to exploit the linearity of 
the equation for the z-component of momentum in LVLH coordinates, namely 

4 = a (2.56) 

Dropping this equation from (2.54), one obtains a reduced system of equation with 
Kronecker indices = {3, 3, 2, 1, . . . , 1} that is linearizable in conjunction with 
(2.56). The reduced system of equations is linearizable with the set of functions 

yu(x) = h\ 
x) = h.% 

V 3 i{x) = 93 

yn(x) = for i G I* 

The remaining state variables are computed by Lie differentiation to obtain the fc th 
time derivative 

y ik {x) = L/yi(jfe_i)(x) = L* -1 y;i(:c) (2.58) 

for i = 1,2,3 and k < K t *. The linear change of coordinates in the input space is 



(2.57) 
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computed from (2.43) by defining 


A{x) 



gl/13 

d$i 

00m 


g&u 

d$i 

00m 

gm 



00m 

/>tji 

• • * Si^m 


••• J 


(2.59) 


where the notation i*j denotes the y th element of the index set I*. A second set of 
Lie derivatives is computed according to 


Wi(x) = L K /yn(x) (2.60) 

for i = 1,2,3 and u\-(x) = 0 for 4 < i < m. The linearizing feedback and input 
coordinate change are given by 

a(x) = — A~ l (x)w(x) and /3(x) = A'" l (x) (2.61) 

The state variables of the linearized system consist of the components of mo- 
mentum in the x and y directions in LVLH coordinates, the first and second time 
derivatives of these momentum components, the momentum in the z direction in 
LVLH, the third component of the attitude quaternion, q 3 , its time derivative, g 3 , 
and a subset of the available gimbal angles. The control variables are the third time 
derivatives of the momentum components /if and /if, the second time derivative 
of q 3 , and a subset of the original gimbal rates. In terms of the notation from the 
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previous section, (z,v) are given by 


A 
Z = 


h f 
h* 

4 
4 

9 3 
93 

\ r / 


and v = 


sf' 

i if 


93 

\ ** / 


( 2 . 62 ) 
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Chapter 3 

Application of Feedback Linearization to the 
Steering Problem 


In the introduction to this thesis, CMGs were described as actuators capable of 
applying torques to a spacecraft for the purpose of attitude control. From the 
derivation of the equations of motion in Chapter 2, it can be seen that a system 
of CMGs is effectively a momentum storage system. By manipulating the rotors 
of individual gyroscopes - through the gimbal drives - it is possible to change the 
rotation rate of a spacecraft by changing the net momentum of its CMG system. 
Although other approaches exist [38] [4], spacecraft autopilots and control systems 
usually command their actuators to apply a desired torque, T det , to the vehicle to 
regulate its attitude. Using CMGs to generate the torques requires solving (1.2) for 
a set of gimbal rates which satisfies 

r c = r de , - (ui B x h cmg (6)j = T>[h cmg (9)]e (3.1) 

The problem of finding a set of gimbal rates is complicated by the fact that the 
number of available gimbals, m, is greater than the number required to provide 
three-axis control of the spacecraft. Also, there are singular configurations of the 
rotors for which the rank of the matrix D[h CTng ] is less than three, implying that 
it is not possible to apply torques to the vehicle in all directions. The problem of 
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finding a set of gimbal rates 9 which satisfies (3.1) and moves the rotors away from 
singular configurations is called the CMG steering problem. 

The CMG steering problem is related to the inverse kinematic problem en- 
countered in robotics. The inverse kinematic problem for a robot involves finding 
trajectories that take an end-effector from an initial position and orientation to 
another specified configuration. The related CMG steering problem is to move in- 
dividual rotors to obtain a change in total stored momentum. The relationships 
between the momentum and gimbal angles of a CMG system, (h,0), and the end- 
effector position and joint angles of a robot, (x,g), are shown in Table 3.1 [2]. The 
robotics literature in particular contains a variety of approaches to this problem. 
Most approaches used in solving the CMG steering problem use the Moore-Penrose 
pseudoinverse solution of (3.1) to compute a set of gimbal rates that produces the 
desired torque, and then add a set of null rates. The null rates are chosen to re- 
distribute the CMG rotors toward more desirable configurations while producing 
no torque. This approach is related to the resolved motion method encounted in 
the robotics literature [26]. In both the CMG and robot problems, a performance 
objective or optimality criterion is used to characterize the desirability of a given 
state or control. In the CMG problem, a set of gimbal rates is chosen by optimizing 
the projection of the rates onto the gradient of the objective function subject to 
the constraint that a specified torque is produced. This process satisfies the torque 
requirements while acting to avoid singular configurations. 

In this chapter, an approach to CMG steering is developed that employs the 
linearizing transformation developed in Chapter 2 for the spacecraft without ex- 
ternal torques. Although the system to be controlled has been made linear by the 
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MANIPULATOR 


SG CMG SYSTEM 


<£> 

Position x = x (q) Momentum h_ 3 h{6) 

Velocity = J(q) q Torque h = J(0)£ 

Acceleration x ■ J (cj) cj + J (q) q Torque h s J(0)0+j{0)£ 

Singularity 

No motion possible in a certain No torque possible In a certain 

direction direction 


Table 3.1: Kinematical Quantities in CMG and Robotic Problems. 


57 



transformation (2.52), the steering problem has not been solved. The first difficulty 
in implementing a steering law is due to a property of the transformation itself: The 
attitude dynamics of the spacecraft are made linear by using three CMG gimbals to 
apply torques as needed to eliminate cross-axis coupling of the attitude dynamics 
and to decouple the spacecraft from the remaining gimbal drives. In effect, the re- 
dundant actuators have been decoupled from the spacecraft. The second difficulty 
is that there are no provisions for steering the CMGs away from undesirable config- 
urations. The linearizing transformation solves a kinematical problem; it does not 
help with the qualitative requirements of the CMG steering problem. If it were not 
for the kinematic redundancy in the CMG system, the solution to the linearization 
problem would be adequate and unique; no steering problem would exist. This 
feedback approach has already been applied to a number of robotic problems which 
are not kinematically redundant [39], [24]. 

Note that the feedback part of the linearizing transformation (2.52) yields a 
parameterization of all sets of gimbal rates that produce a specified attitude ac- 
celeration. Once the first three elements of the control vector of the linear system 
have been fixed, the remaining elements can be selected arbitrarily without chang- 
ing the commanded acceleration. The remaining elements are used to determine 
null motions. The parameterization effectively eliminates the constraint associated 
with the steering problem and reduces the number of variables that need to be con- 
sidered in an optimization three. This parameterization can be applied to existing 
steering algorithms to simplify the numerical optimization problem and to reduce 
the number of required calculations. 

A number of algorithms for CMG steering and robot control select a solution 
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by using a generalized inverse of the form 


9 = C'r e + (J m - C*C)w (3.2) 

where C = D^] and C* = C T (CC T )~ 1 is the Moore-Penrose pseudoinverse of 
C. This expression is similar to the parameterization obtained from the linearizing 
feedback in that the second term in (3.2) projects an arbitrary vector w G R m 
into the null space of C . Because w is an m-vector, the expression (3.2) does not 
give a minimal characterization of the null motions; a fixed w does not map to a 
unique set of gimbal rates. In addition, the number of arithmetic computations 
required to evaluate (3.2) is greater than that required to evaluate the linearizing 
transformation. The number of computations involved in calculating the matrix A 
in the linearizing transformation (2.52) and the matrix C are essentially the same. 
Inversion of the matrix A can can be reduced by column operations to inversion 
of a 3 x 3 matrix and multiplication of the inverse by a 3 x (m — 3) matrix; more 
computations are required just in the computation of (7*. 

In the next section, an approach to distributing control effort among a group of 
redundant actuators is presented. This approach uses the linearizing feedback to 
state the distribution problem as a normed approximation problem. The following 
section then develops this approach for solving the CMG steering problem. This 
strategy has a conceptual advantage over conventional approaches to CMG steering 
in that developing a steering law requires developing an algorithm for computing 
desired gimbal rates, rather than the standard method of developing a function 
whose gradient yields the desired rates. In addition to being more direct, such an 
algorithm can be programmed to recognize problem situations and take specific 
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corrective measures. Finally, several examples are presented which show the per- 
formance of a steering law developed for the parallel mounted CMG configuration 
proposed for the NASA space station. 


3.1 Formulation of the Normed Approximation Problem 

In the proof of the feedback linearization theorem presented in Chapter 3, it was 
assumed that the dimension of the distribution Q 0 was equal to m, the number of 
inputs in (B.52). If there are redundant inputs, it is possible to find an invertible 
matrix function H(x) such that 


G{x)H{x) 


9l(») ••• Sm( X ) 0 ■■■ 0 


(3.3) 


Under this coordinate change, it is only necessary to find a transformation that 
linearizes the system with inputs ft. If the input transformation (a, j3) solves the 
reduced problem, then the feedbacks given by 



a 


0 o 

a = H(x) 


and /? = H( x) 


0 


0 Im—TTi 


solve the original linearization problem. 

Clearly, the columns of the matrix H(x) are related to the right eigenvectors 
of the matrix G(x). The solution to the linearization problem given by (3.3) may 
be of little practical value because the redundancy has been eliminated by the 
transformation H. The last m — rh control variables have their influence in the 
null space of G(x) - they do not influence the nonlinear system being controlled. 
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It is quite possible that the linearization is such that certain actuators are used 
excessively for control purposes, while other actuators are not used at all. In effect, 
the burden of control is placed on only a few of the available actuators, and these 
actuators may not have the physical authority to obtain the desired results. As a 
result, linear control laws based solely on the linear system may not yield reasonable 
performance. In fact, a similar situation occurs with the two spacecraft models that 
have been linearized. 

Redundancy appears in the spacecraft problems in a slightly different manner 
than has just been described. Even though only three CMG gimbals are required for 
control of the spacecraft, they appear in the linearized coordinates as independent 
inputs to the nonlinear system. In other words, the distribution Q 0 has maximal 
rank. The reason that the actuators are independent is that there are dynamics 
associated with the gimbal drives; however, the difficulties that have been described 
for systems with redundant inputs still occur. The feedback transformation uses 
three of the gimbals to control the spacecraft, and the remaining gimbal rates can be 
commanded arbitrarily. Generally, three of the gimbal drives will not have sufficient 
authority to perform attitude control functions alone. This is especially true if a 
number of other CMGs are being freely gimballed, or if the three gimbals are in 
a poor configuration for producing acceleration in the desired direction. Before 
control laws based on the linearizing transformations can be effectively applied, it 
is necessary to develop an approach that distributes the control effort required to 
produce a specified attitude acceleration among the available actuators. 

Redundancy arises in the spacecraft problem because there is a subsystem in the 
model equations that is to be controlled. Here it is if interest to control the space- 
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craft; the fact that there are dynamics associated with the actuators is incidental in 
this case. For the purpose of control, it is desirable to have the spacecraft attitude 
dynamics appear independently of other system dynamics. More specifically, the 
linearization problem could have been stated with the additional requirement that 
yn(x) = qi for i = 1, 2, 3 in (B.89). This requirement necessarily implies that = q^ 
for i = 1,2,3. The linearization problem with this type of restriction on some of 
the state variables has been solved in [24]. The linearizing feedback for this class of 
problems has the form 


u 




02 


V 2 


(3.5) 


where the set of control variables (ui,/?i) influence the subsystem whose dynamics 
are of particular interest to control. The remaining pair (v 2 ,/3 2 ) can be associated 
with other subsystems, or with the dynamics of actuators that are not required for 
controlling the important subsystem. 

A natural way to resolve the redundancy issue would be to choose a set of inputs 
that requires a minimal effort to produce a desired influence on the controlled sub- 
system. This would be adequate in the case where there are no dynamics associated 
with the actuators. In more general situations, however, it may be desirable to drive 
the state of the actuators to a particular value. Rather than simply minimizing the 
control effort, there may be a desired value Ud that the input should take. The 
problem of computing an input variable u that is close to the desired value can be 
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stated as a normed approximation problem 

min Hu,* — u|| (3.6) 

u 

subject to the constraint that the controlled subsystem is influenced in a speci- 
fied way. Note that (3.5) provides a way to eliminate the constraint in the above 
problem. The variables Vi were chosen because they are time- derivatives of system 
parameters which are to be controlled. Choosing V\ appropriately yields a set of 
inputs to the nonlinear system that produces the desired result. 

Letting a = uj — a — /3iVi, the contrained optimization problem can be restated 
as an unconstrained optimization in terms of the variable v 2 


min ||o — f3 2 v 3 \\ (3.7) 

V2 


From a practical point of view, it is useful to impose limits on the individual inputs 
of the form 


— u 


peak 


< U{ < U 


peak 


(3.8) 


where u l peak is a physical limitation associated with device i. Fast numerical methods 
exist for solving this problem in three important norms: The 1-norm defined by 


Nil = it\ X i\ 

1=1 

the 2-norm is the usual Euclidean norm defined by 



(3.9) 


(3.10) 
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and the oo-norm is defined to be 


ML = max l x *'l ( 3>11 ) 

Xi 

In the following three subsections, the 1-norm and oo-norm problems are formulated 
as linear programming problems, and the 2-norm problem is stated as a quadratic 
programming problem (numerical methods are discussed later). 


3.1.1 Solution of the 1-norm Approximation Problem 

To formulate the H*^ optimization problem, a new set of variables are introduced 
by defining the vector y as 


y = u d - a - favi - /3 2 v 2 (3.12) 

The new variables yi represent the error in approximating the desired rates and are 
written in the form 

Vi = vt - Vi (3.13) 

with yt,yr ^ 0 f° r 1 < i < m. The linear program that solves the optimization 
problem is 

m 

min ||ud-u|| x = min J^yt + yf (3-14) 

v 2 v 2 ,y + ,y~ i = i 
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subject to the three sets of m constraints 


P2V2 + y + - y“ = « 

a + P 1 V 1 + P 2 V 2 < u peak ( 3 . 15 ) 

-a-fJiVi-p 2 V2 < Upeak 

In this formulation of the problem, the decision variables yf , are non-negative, 
and the (u 2 )« are free. 


3.1.2 Solution of the 2-norm Approximation Problem 


The minimization in ||-|| 2 can be formulated as a quadratic program in a straight 
forward manner by defining 


min Hud — uHj 

v 2 


min {ud — u) T (ud — u) 

v 2 

min (a T a — 2a T f3 2 v 2 + fl 2 v 2 S ) 
v 2 


( 3 . 16 ) 


subject to the two sets of m constraints 


a + /3ivi + f3 2 v 2 < u peak 
-a - (3 1 v l - f3 2 v 2 < u peak 


( 3 . 17 ) 


In this formulation, the decision variables can be positive or negative, and the 
constraints impose limits on the allowed values. 
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3.1.3 Solution of the oo-norm Approximation Problem 

The ||*|| 00 optimization problem is formulated by first introducing a new variable, 
A, such that 

- A < ^ - (0 2 v 2 )i < A (3.18) 

for 1 < i < m. This is equivalent to 

k-(/?2«2)i| < A Vi (3.19) 

and A is a bound on ||*|| . The optimization problem is to minimize A 

min ||ud — u|| TO = min A (3.20) 

v 2 v 2 

subject to the four sets of m constraints 

A + 0 2 v 2 > a 

A — (3 2 v 2 > —a 

(3.21) 

- 0 2 v 2 > -u peak + a +/3 1 v 1 
(3 2 v 2 > - u peak - a - 0 lUi 

Since revised simplex algorithms require more computational overhead to deal with 
a given number of constraints than a given number of decision variables, it is more 
efficient to solve the dual problem associated with (3.20,3.21). The dual linear 
program is stated in terms of the variable w = [uq, . . . ,^4 m ] T a s 

T71 

max J2 aiVJ i ~ a < w 2 i + {-u peak + a + fav^iWzi - (u pea fe + a + 0\V X )iW Ai (3.22) 

W i=l 
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subject to the m + 1 constraints 


2m 

XX ^ 1 

i= 1 

'■ m 

XI [&1.M ( w i _ w 2i + ™3j - «>4i) = 0 for 1 < i < m 

i = i 


( 3 . 23 ) 


The problem is formulated such that all decision variables are non-negative, i.e. 

> 0 V i. 


3.2 A Steering Algorithm Based on Normed Approxima- 
tion 

In this section, an algorithm for computing a desired set of gimbal rates for a group of 
double-gimballed CMGs is proposed. The algorithm is based on computing a desired 
rotation vector for each CMG. The direction and magnitude of the rotation vector 
of a CMG rotor correspond to the axis and rate of rotation, respecively, of the rotor 
due to gimbal motion. The rotation vectors are chosen to move the CMG rotors 
to maximize the amount of momentum which can be projected instantaneously in 
any direction, to avoid line-up with the commanded torque vector, and to keep 
inner-gimbal angles small. The desired rotation vector for a CMG is computed as 
the sum of a set of rotation vectors each of which individually attempts to achieve 
one of these objectives. From the desired rotation vector, desired gimbal rates are 
derived for each CMG. This section concludes with a brief explanation of how the 
proposed steering law is incorporated into an attitude control system. 
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3.2.1 Computation of a Desired Set of Gimbal Rates 


One method for maintaining good three-axis controllability with double-gimballed 
CMGs is to keep the angles between rotors as large as possible. Achieving this goal 
corresponds to maximizing the minimum torque that can be projected in any direc- 
tion and corresponds to keeping three-axis controllability of the spacecraft high in 
all directions (this is not equivalent to gain maximization [7]). Since singular states 
of double-gimballed CMG systems are always associated with a mutual alignment 
of rotors, this is also equivalent to avoiding singularities. For each CMG, the term 
which achieves this in the rotation vector computation is 


<r\ = — [cos 1 (h t • h*) — 7rl unit{h x x W) (3.24) 

i*i * 


where u p denotes the peak gimbal rate, h l denotes a unit vector along rotor axis 
z, and the cross product is normalized to have unit norm. Rotating each CMG in 
this way will move the rotors to a configuration in which the rotors are distributed 
symmetrically about the total CMG momentum vector. The desirability of such 
configurations has been pointed out in [4]. 

If a rotor is allowed to align with the commanded torque, then that rotor cannot 
contribute to the torque. The term which steers rotors away from such configura- 
tions is 


if cos x (h l • f ) > k T 


^2 = 1 


(3.25) 


-jZ [cos 1 (h l • f) — fc r j unit{h l x f) if cos 1 (h t ■ f) < k T 
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where f denotes a unit vector in the direction of the requested torque, k T is a 
constant, and the cross product is normalized as before. Some investigators have 
suggested steering toward the negative of the desired torque to avoid singularities [8]. 
To avoid numerical difficulties, an arbitrary gimbal axis is used in the neighborhood 
of an exact line-up of vectors in (3.24) and (3.25). 

To simplify notation in defining the next two contributors to the rotation vector, 
the function $ is defined by 

r 

0 if |t|^ c 

i(a,b,c,d) = | jjfii if c < ||| < d (3.26) 

1 ^ l > d 

where (a, 6, c,d) are constants with \a\ < b , and 0 < c < d < b. 

It is desirable to keep the inner gimbal angles small, both to avoid potential 
hardware stops, and because the torque due to the outer gimbal rate is proportional 
to the cosine of the inner gimbal angle measured from the reference angle. A rotation 
vector which reduces the inner gimbal displacements is 

4 = -t (3.27) 

where Q{ denotes the displacement of gimbal t, 0 Btop is the displacement to a gimbal 
stop, fcj nner , k l ™ er are constants, and <Tj n ncr is a un it vector along the inner gimbal 
axis. 

In studying several maneuvers, it was found necessary to include a term in the 
calculation of the desired rotation vector which would rotate the outer gimbal angle 
if the inner gimbal was being forced against a stop. Associated with the inner 
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3.2.2 Implementation of the Proposed Algorithm 


In this section, an algorithm is described which implements a CMG steering law 
based on the normed approximation approach to redundancy management. The 
program has been implemented in FORTRAN on a VAX computer along with a 
rigid body simulation of a spacecraft. A flow chart of the steering algorithm is shown 
in Figure 3.1. The inputs are the current CMG gimbal angles, the spacecraft atti- 
tude quaternion, the derivative of the quaternion, and the second derivative of the 
quaternion, which is determined by an autopilot or control law. Since quaternions 
are algebraically related to any set of attitude variables, the proposed implementa- 
tion is quite general. 

The steering algorithm first tests the CMG gimbal angles to see if the CMG con- 
figuration is near a singular state. If the gain function is above a certain threshold, a 
routine that monitors the transformation is called. The routine MONITOR checks 
that the three gimbal rates that do not appear in the linearized control vector are 
capable of producing torques in all directions (these three gimbal rates are used to 
implement the linearizing feedback). Their ability to control the spacecraft is deter- 
mined by the determinant of the matrix formed by the three columns of the matrix 
Dfhcmg] corresponding to these gimbal rates. When this determinant is small, the 
matrix that must be inverted during the computation of the linearizing transfor- 
mation is nearly singular. If the determinant is below a certain threshold, a new 
combination of three gimbal angles is found which produces the largest determinant, 
and these three gimbal rates will not appear in the transformed control. 

If the CMG rotors are near a singular configuration, then there will be no com- 
bination of three CMGs that can produce torques adequately in all directions. Most 
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gimbal vector of each CMG is the rotation rate 


a\ = -u peak m,0>top,k? Ur ,kr' r ) Pouter (3.28) 

where 0{ and 0, top are as before, k^ uteT , fcf uter are constants, and a , outer is a unit vector 
along the outer gimbal axis. The initial sign of this rotation is chosen to move the 
rotor awav from the inner gimbal rotation plane, and is not changed until <r\ is set 
to zero. The need for such a term will be explained in connection with one of the 
examples. 

For each CMG, the desired rotation vector is found by taking a weighted com- 
bination of the four rotation vectors defined above. 

< = £, k cr} (3.29) 

J = 1 

From this composite rotation vector, it is possible to compute a unique set of gimbal 
rate commands that produce the desired rotor motion for the given CMG by solving 
the linear equation 

4 x V = uf' r (dinner X V) + U?'" (<T out „ X V) (3.30) 

for the desired inner and outer gimbal rates. Proceeding in this way for each CMG, 
a set of desired gimbal rates, tij, is computed, and the normed approximation 
strategies of Section 4.1 can be used to calculate actual gimbal rates based on these 
desired rates. In the following subsection, an implementation of a steering law based 
on this desired rate computation and the approximation algorithms is described. 
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Figure 3.1: Flow Chart of Attitude Control Algorithm 
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steering algorithms do not consider this situation at all; even though no steering 
algorithm can really guarantee absolute singularity avoidance. The most common 
symptom of failure is termed gimbal lock , because the rates required to produce a 
desired torque become so large that no null motion is achieved. For CMG systems 
consisting of single-gimballed devices, there are singular states for which there is no 
possible escape through null motion, and the CMG system will lock. It is neces- 
sary that a steering algorithm continue to operate in a sensible manner under these 
circumstances. 

Near singular states, the linearization problem is not solvable as stated, and 
the problem must be modified so that the algorithm can proceed. The numerical 
difficulty that arises in computing the transformations is due to the fact that the 
matrix A{ x) in (2.52) is nearly singular. By exchanging columns, this matrix can 
be brought into the form 

Aj(x) A 2 {x) 

0 An- 3 

where A\ is a 3 x 3 matrix. The MONITOR routine attempts to pick a set of 
three gimbal angles that guarantee numerical invertibility of Ax(®); however, near 
a singular state, no such set of gimbals exists. The matrix Ai can be decomposed 
into 

At = V t AV (3.32) 

where A is the diagonal matrix of eigenvalues of A l5 and V is the matrix whose 
columns are the orthonormalized eigenvectors. The problem of computing the lin- 
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earizing feedback, (2.52) can be transformed to the form 


V 

0 


Vq 


A V 

v a 2 



V — i 


► = 



0 

Im- 3 


6' 


0 

Im- 3 


{w(z) + u } 


(3.33) 


where the third row in the matrix on the right is nearly zero - assuming the eigen- 
values decrease along the diagonal of A. Eliminating this row from the equation 
and increasing the dimension of the identity matrix by one yields the feedback 
transformation 
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Near a singular state, there is a direction about which the CMGs can produce no 
acceleration. This modification to the algorithm projects the desired acceleration 
onto the space orthogonal to the singular direction, and control of the spacecraft 
about the singular direction is replaced by an extra gimbal rate in the transforma- 
tion. Indeed, the spacecraft will not be controlled in the singular direction; however, 
control will be maintained along the other axes. Controllability will be recovered in 
the singular direction after the singular orientation has been eliminated. Depending 
on the relative orientations of the singular direction and the torque direction, sin- 
gularity may be a problem for only a very short time. This modification guarantees 
that the linearization problem is always solvable. 

After the program resolves potential computational difficulties, the linearizing 
state and feedback transformations are evaluated in the routine TRANSFORMS. 



The routine CMG-DISTRIB is then called to compute the desired set of gimbal 
rates based on the current state of the CMG system. The routine APPROX calls 
the appropriate subroutines for formulating and solving the desired normed approx- 
imation problem. The 1-norm and 2-norm approximation problems are solved using 
the subroutines DLPRS and QPROG, respectively, from release 10.0 of the IMSL 
subroutine library. The oo-norm approximation problem is solved using subroutine 
ZX3LP from release 9.2 of the IMSL library. The output of the steering program 
is a set of gimbal rates, u which produce the commanded q. 

This steering algorithm has been implemented as part of a simple simulation of a 
rigid spacecraft. The equations of motion are integrated using a fourth-order Runge- 
Kutta integration scheme. Before computing an integration step, a subroutine is 
called that can implement a variety of attitude control functions, and then the 
CMG steering algorithm is called. The chosen gimbal rates are then passed to a 
subroutine that executes a single integration step. The following section presents 
results from simulations using the program described in this subsection. 

3.3 Simulation Examples 

In this section, several simulations are presented to demonstrate the behavior of 
the proposed steering algorithm performing a variety of attitude control functions. 
The first three examples illustrate the characteristics of the gimbal angle trajectories 
generated by the three normed approximation algorithms. It will be shown that the 
1-norm and oo-norm algorithms occasionally exhibit undesirable oscillations in the 
present application. The 2-norm algorithm generally yields well behaved solutions 
and is used in the majority of the examples. The fourth example demonstrates 
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the need for the fourth term in the rotation axis computation, (3.28). This is an 
example of a situation in which the steering algorithm must deal explicitly with 
a problem that conventional gradient based methods would fail to recognize. The 
next two examples show how the algorithm redistributes the CMG rotors when 
they are initialized in poor configurations. In one of the examples, the algorithm is 
initialized in a singular configuration. The final example demonstrates the response 
of the vehicle/CMG system when a large onboard disturbance is applied. The 
disturbance is large enough to rate saturate several of the CMGs, and the example 
demonstrates the advantage of including rate limits on the individual gimbals. 

The inertia properties of the vehicle are listed in Table 3.2 and represent the 
rephased dual-keel space station, after assembly sequence flight 7 [40]. Attitude 


Inertia Matrix slug-ft 2 

2.425 x 10 7 

-1.917 x 10 5 

3.060 x 10 5 

-1.917 x 10 5 

3.556 x 10 6 

2.333 x 10 5 

3.060 x 10 5 

2.333 x 10 5 

2.515 x 10 7 


Table 3.2: Inertia Properties for Re-phased Dual-keel Space Station. 

control torques are provided by four parallel-mounted double-gimballed CMGs. In 
the parallel mount configuration, the outer gimbals are aligned with the vehicle pitch 
axis and the inner gimbals are aligned along roll. Each CMG rotor possesses 5000 
ft-lb-aec of angular momentum and can be gimballed about inner and outer mounts 
at peak rates of 5 deg/sec. Inner gimbal displacements are limited to less than 
90 deg and outer gimbals may rotate continuously. Figure 3.2 shows a schematic 
diagram of the vehicle along with the body reference frame, and the CMG mounting 
configuration is shown in Figure 3.3. The set of constants used by the steering law 
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Rotation Vector Constants 

fcj = 0.2165 


k 2 = 1.0 

k r = 1.57 

k 3 = 1.0 


1 fe 4 = 1.0 

osbxhhh 


Table 3.3: Constants Used in Rotation Vector Computation. 

to compute the desired gimal rotation vectors are shown in Table 3.3. Although 
these simulations employ parallel CMG mounting, the structure of the steering law 
by no means requires this, and is ammenable to any mounting configuration. 
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Figure 3.2: Schematic of Rephased Dual-Keel Space Station; Flight 7. 



Figure 3.3: Reference Configuration for Four Parallel Mounted CMGs. 


78 




3.3.1 CMG Saturation Along Pitch 


The first set of simulations illustrate the characteristics of the gimbal angle trajec- 
tories generated by the three normed approximation algorithms. The CMG system 
is initialized in a zero momentum state with gimbal angles (in degrees): 


1 31 

^inner 


5.73 


= o.o 


£2 

''inner 

03 

^inncr 


-5.73 «L» = -90 

and 

-5.73 el ultr = 90 


9\ nn „ = 5.73 d* outer = 180 


(3.35) 


(Note that 0.1 rad. = 5.73 deg.). In each of the three cases in this subsection, 
the CMGs are commanded to produce a torque along the negative pitch axis of the 
vehicle. This causes the net CMG momentum to saturate along the positive pitch 
axis. Figure 3.4 shows the rotor trajectories projected onto a unit sphere for the 
case in which the 2-norm approximation is used to compute the gimbal rates, and 
the torque produced is a constant -1 95 ft-lb. In the starting configuration, the four 
momentum vectors (CMG rotors) nearly align with the positive and negative roll 
and yaw axes. The requested torque causes the rotors to project momentum along 
the positive pitch axis and the rotors align along this direction. The gimbal angles 
and gimbal rates are shown in Figure 3.5 and Figure 3.6, respectively. For most of 
the simulation, the desired inner gimbal rates are positive since saturating along the 
pitch axis requires pushing the inner gimbals against their (negative) stops. Since 
inner gimbals 2 and 3 have a negative dispacement initially, their commanded rates 
are smaller until all inner gimbal angles are approximately the same. After about 
85 sec ., the outer gimbals are forced to rotate by the fourth term in the rotation 
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axis computation (3.28). This causes the rotors to saturate unevenly during the 
last 10 sec . 

In Figure 3.7 and Figure 3.8 the gimbal angles and rates are shown for the 
simulation in which the 1-norm approximation algorithm is used. Again, the com- 
manded torque is - 195 ft-lb . At the beginning of the simulation, the inner gimbal 
angles converge linearly until the angles are nearly the same. The inner gimbal 
angles then oscillate about a nominal trajectory. For a time, Figure 3.8 indicates 
that these oscillations are caused by the inner gimbal rates being switched between 
the positive and negative peak rates. This undesirable behavior occurs when the 
gimbal angles are near an optimal configuration, and is due to the fact that the 
approximation algorithm is computing an instantaneously optimal solution, not a 
continuous solution. A CMG configuration is referred to as optimal when the de- 
sired gimbal rates have almost no null component. In other words, redistribution at 
constant momentum can not improve the orientation, hence the only component in 
the calculation of the desired rates reflects a desire to transform momentum in the 
opposite direction required by the commanded torque. When the oscillations begin, 
one pair of gimbals is displaced slightly more toward the inner gimbal stop than 
the others, and hence, has a larger positive desired gimbal rate reflecting a desired 
to reduce the inner gimbal angle. The optimal solution attempts to approximate 
these rates accurately and allow the rates of the other pair of gimbals to be large 
and negative. During the next computational cycle, the situation is reversed, and 
the oscillations result. 

A number of attempts were made to eliminate the oscillations that occur when 
the 1-norm approximation is used. The most successful used a first-order linear 
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filter to to average the gimbal rates computed by the approximation algorithm to 
produce a smoother response. The filter suppresses the oscillations, but imposes an 
unacceptably low bandwidth when torque levels are high. 

In Figure 3.9 and Figure 3.10, the gimbal angles and rates are shown for the 
simulation in which the co-norm algorithm is used. In order to avoid causing a 
gimbal to be pushed against a stop before the end of the simulation, the requested 
torque is -1 90 ft-lb. Like the previous two cases, the inner gimbals converge on the 
same value as they begin to saturate. The co-norm alogorithm allows significant 
motion of the outer gimbals, even though this motion cannot contribute to the 
torque applied to the spacecraft. The excessive outer gimbal motion occurs because 
the optimization is attempting to minimize the largest absolute error in approxi- 
mating the desired gimbal rates. Nominally, the desired rates for the outer gimbals 
are zero, while the desired inner gimbal rates increase as the simulation proceeds. 
Since the error in approximating the desired inner gimbal rates becomes large, it is 
clear from the definition of H'H^ that the outer gimbals can be commanded to rotate 
without influencing the optimal value of the approximation. In fact, if it is possible 
to reduce the peak error by any amount using the outer gimbals, perhaps due to 
numerical round-off, they will be employed at the largest possible rate consistent 
with the norm of the optimal solution. Evidence of this can be seen by noting that 
the peak values of the outer gimbal rates approximately correspond to the value of 
the desired rate for the inner gimbal angles shown in Figure 3.11. 

Both the 1-norm and oo-norm approximation algorithms can exhibit undesirable 
oscillations under certain circumstances, and it is not desirable to command CMG 
gimbals to follow these kinds of rates. The chattering generally occurs around a 
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near-optimal configuration, this undesirable behavior does not occur when there is 
an unambiguous policy to follow in redistributing the rotors. For example, Fig- 
ure 3.7 shows the initial convergence to an optimal configuration is quite fast using 
the 1-norm approximation, indicating that both the 1-norm and oo-norm approx- 
imations may have useful applications in certain circumstances. For example, the 
1-norm is particularly relevant when the control is related to a displacement, and 
the oo-norm solution is useful when it is desired to minimize the maximum effort 
required from each individual actuator. Because of its intrinsic smoothness, the re- 
maining simulations were all performed using only the ||*|| 2 optimal approximation. 
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ROLL 



PITCH 


Figure 3.4: Saturation of CMG Rotors along Pitch Axis (668). 
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INNER-G IMBfiL DISPLACEMENTS 

o 



OUTER-G IMBfiL DISPLACEMENTS 

o 



Figure 3.5: Gimbal Angles; Momentum Saturation along Pitch, 2-norm (668). 
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DEGREES / SECOND DEGREES / SECOND 

- 7.5 - 5.0 - 2.5 0.0 2.5 5.0 7.5 - 1.0 - 3.0 - 2.0 -|.0 0.0 


INNER-GIMBRL RATES 



OUTER-GIMBAL RATES 



Figure 3.6: Gimbal Rates; Momentum Saturation along Pitch, 2-norm (668). 
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INNER-GIMBRL DISPLACEMENTS 



0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 100.0 110.0 


TIME (SEC0N0S) 



Figure 3.7: Gimbal Angles; Momentum Saturation along Pitch, 1-norm (669). 
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DEGREES / SECOND DEGREES / SECOND 

- 2.0 0.0 2.0 4.0 6.0 - 7.5 - 5.0 - 2.5 0.0 2.5 5.0 


INNER-GIMBflL RRTES 
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0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 100. 0 110.0 

TIME (SECONDS) 


Figure 3.8: Gimbal Rates; Momentum Saturation along Pitch, 1-norm (669). 





INNER-GINBRL DISPLACEMENTS 



0UTER-GIM8RL DISPLACEMENTS 



Figure 3.9: Gimbal Angles; Momentum Saturation along Pitch, oo-norm (731). 
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DEGREES / SECOND DEGREES / SECOND 
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TIME (SEC0N0S) 


0UTER-9IMBRL RATES 



0.0 10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 90.0 100.0 110.0 


TIME (SECONOS) 


Figure 3.10: Gimbal Rates; Momentum Saturation along Pitch, oo-norm (731). 





DESIRED INNER-GIM8RL RRTES 



DESIRED OUTER-GIMBRL RRTES 



Figure 3.11: Desired Gimbal Rates; Momentum Saturation along Pitch, oo-norm 
(731). 
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3.3.2 CMG Saturation Along Roll 


In the second example, a torque of -1 00 ft-lb is commanded along the negative roll 
axis of the vehicle, and the initial configuration is the same as in the previous case. 
This torque will cause the CMG rotors to align along the positive roll axis. Because 
of the particular symmetry of the situation, the desired rotation vector of rotor 2 
consistantly lies along its inner gimbal axis. As a consequence, this gimbal will 
eventually be rotated against its stop. The outer gimbal rotation term (3.28) was 
added to the rotation vector computation to prevent this from happening. 

Figure 3.12 and Figure 3.13 show the gimbal displacements and gimbal rates 
from the simulation. At about 125 sec. into the simulation, the outer gimbal of 
rotor 2 begins to rotate at its peak rate. Initially, the inner gimbal axis is aligned 
with the negative yaw axis of the vehicle. As the outer gimbal angle rotates, the 
inner gimbal axis eventually has a component along the positive yaw axis and the 
gimbal rate must change sign in order to continue projecting more monentum in the 
positive roll direction. This causes the sudden change of sign in the inner gimbal 
rate. Figure 3.14 shows that rotating the outer gimbal causes a drop in the gain 
function shortly before t=125 secs. Steering laws based on gain optimization would 
never choose to rotate the outer gimbal in this situation, and the inner gimbal 
would be forced against a stop. A similar statement applies to steering laws that 
attempt to avoid rotor line-up. This is an example of a situation in which the 
steering algorithm must deal explicitly with a problem situation. Note also that the 
approximation algorithm is bounding the rate of outer gimbal 3 at its peak value. 
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INNER-GIMBAL DISPLACEMENTS 



OUTER-GIMBflL DISPLACEMENTS 



Figure 3.12: Gimbal Angles; Momentum Saturation along Pitch, 2-norm (674). 
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Figure 3.13: Gimbal Rates; Momentum Saturation along Pitch, 2-norm (674). 
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Figure 3.14: Gain Function during Momentum Saturation along Pitch, 2-norm 
(674). 
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3.3.3 CMG Redistribution from Poor Configurations 

In the next two simulations, an angular rate change is commanded with the CMGs 
beginning in an initially poor configuration. The attitude control law commands an 
angular acceleration proportional to the error in the angular rate. The vehicle rate 
response to a step change in rate of -0.23 deg /sec in pitch is shown in Figure 3.15. 
The rotors are initialized in a zero momentum state with rotors 1 and 2 aligned 
with the positive roll axis, and rotors 3 and 4 nearly aligned with the negative roll 
axis; this configuration is nearly singular. The gimbal angle responses are shown in 
Figure 3.16. The inner gimbals must project momentum along the positive pitch 
axis to achieve the rate change, but the outer gimbals move until they are evenly 
distributed about the net CMG momentum vector in order to alleviate the singular 
condition and attain a best orientation. In Figure 3.17 the gain function is shown. 
The gain is initially almost zero and increases rapidly to about 80% of its maximum 
value (recall that gain is not being optimized). It then decreases as the CMGs must 
project momentum in the pitch direction to maintain the rate change. This example 
illustrates that the proposed algorithm functions in the vicinity of a singular state. 
Even if the requested torque had been along the singular direction, null motions 
would drive the rotors to a configuration which could produce the desired torque. 

The same rate response was simulated with the CMGs beginning in another 
poor configuration with net momentum projected in the roll direction. In this 
case the inner gimbal angles are initially zero and the outer gimbal angles are 
0 o = [— ~, — — ~] ra d (rotors 1,2 and 3,4 are initially aligned). The gimbal 

angle histories are shown in Figure 3.18. Once again, the inner gimbals displace 
negatively to project momentum in the pitch direction, while the outer gimbals 
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simultaneously displace to avoid line-up. In this case, the final gimbal rate is not 
zero, since the gimbals must rotate at the same rate and in an opposite sense to 
the spacecraft. This insures that torques are not applied and that the net CMG 
momentum remains inertially constant once the vehicle has achieved the desired 
rate. The signs of the gimbal and spacecraft rates are the same because the outer 
gimbal axis points along the negative vehicle pitch axis. 

A difficulty encountered with algorithms that seek local optima of functions is 
that the algorithm can get stuck in a configuration that is locally optimal. For 
the parallel mount configuration of four CMGs, the technique demonstrated here a 
unique, stable, optimal configuration of the rotors provided that the CMGs project 
momentum in either the pitch or the yaw direction (the placement of the rotors is 
unique, not the location of the rotors by index). For each value of net momentum, 
there is also a stationary point of the steering algorithm that is unstable. If the 
CMG project momentum only in the pitch direction, any optimal configuration can 
be rotated as a whole about the outer gimbal axis and still be optimal. 
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Figure 3.15: Attitude Rate Response to Step Change in Reference Rate, 2-norm 

( 686 ). 


97 



degrees oegrees 

-I SOt 000.0 -SO. 0 0.0 50.0 too. 0 ISO. 0 -50.0-40.0 -30.0 -20.0 -|0.0 0.0 


INNER-GIMBRL DISPLACEMENTS 



OUTER-GIMBflL DISPLACEMENTS 



3 

4 



0.0 10. 0 20.0 30.0 


40.0 50.0 60.0 70.0 

TIME (SECONDS) 


80.0 90.0 100.0 110.0 


Figure 3.16: Gimbal Angles; Rate Change along Pitch, 2-norm (686). 
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Figure 3.17: Gain Function during Momentum Saturation along Pitch, 2-norm 

( 686 ). 
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Figure 3.18: Gimbal Angles, 1 Rate Change along Pitch, 2-norm (688). 
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3.3.4 Attitude Response to Torque Pulses 

In this simulation, a feedback control law is used to regulate the spacecraft at 
a reference attitude, q = (0,0,0). Quaternion and quaternion rate feedback are 
chosen so that the quaternion dynamics are characterized by decoupled, second 
order systems with natural frequencies of 0.1 rad/sec and damping ratios of 0.7. At 
the beginning of the simulation, the torque 


T = \ 


— 700 
— 1100 
-700 


ft -lb 


(3.36) 


is applied to the vehicle for 7.5 sec.. Twenty seconds later, the torque is applied 
in the opposite direction. Momentum conserving impulses of this type are repre- 
sentative of onboard disturbances such as crew motion. The attitude response to 
these torque disturbances is shown in Figure 3.19 and the gimbal displacements and 
rates are shown in Figure 3.20 and Figure 3.21, respectively. During the interval 
between the two torque pulses, the gimbals must remain displaced form the start- 
ing configuration, since the initial disturbance causes a change in the total inertial 
momentum of the vehicle. By the end of the simulation, the rotors return to an 
orthogonal configuration with the outer gimbal angles slightly rotated with respect 
to the starting configuration. Note that compensation of the disturbance torques 
saturated some of the gimbal rates. The approximation algorithm recognizes the 
rate constraints, and uses other gimbals accordingly to prevent their violation. The 
only other steering algorithm to consider rate constraints explicitly is [4]. 
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Figure 3.19: Attitude Response to Opposite Torque Pulses (721). 
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Figure 3.20: Gimbal Angles; Opposite Torque Pulses, 2-norm (721). 
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Figure 3.21: Gimbal Rates; Opposite Torque Pulses, 2-norm (721). 
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Chapter 4 

Application of Feedback Linearization to the 
Momentum Management Problem 


A spacecraft orbiting the Earth is subject to a variety of disturbance torques that 
influence its attitude. In general, disturbances due to motion occuring onboard 
such spacecraft are momentum conserving. For example, torques caused by crew 
motion and rotating machinery do not result in a net change in the total momentum 
of the vehicle. Even though there will be a short term fluctuation in spacecraft 
momentum during the interval between push-off and impact of a crew member, 
the total system momentum remains constant. There are, however, a number of 
environmental torques that do cause a net change in momentum. Small torques due 
to gravity gradient and aerodynamic effects can lead to substantial changes in the 
total momentum of the spacecraft over a long period of time. In certain on-orbit 
applications, large transfers of momentum may occur due to disturbances such as 
docking with other spacecraft or venting of waste gasses. It is therefore necessary 
to examine the long-term momentum requirements of CMG equipped spacecraft to 
insure that the momentum capacity of the CMG system is not exceeded. 

Most approaches to momentum management involve manipulating the attitude 
dependent gravity gradient torque to produce a desired transfer of momentum 
[14], [15]. Some authors have studied the use of appendages to manipulate aero- 
dynamic torques [41], [1]. In general, the dynamics of the spacecraft are linearized 
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or approximated in some way, and a dynamical representation of the CMGs is usu- 
ally not included in the vehicle model. In the following section, the momentum 
management problem is formulated as an optimization problem with an integral 
quadratic perfomance index. The performance index is to be optimized subject 
to the constraint that the nonlinear equations describing the spacecraft attitude 
dynamics (2.35) are satisfied. In (2.35) the control variables are the individual 
gimbal rates, and the momentum management problem formulated here includes a 
dynamical model of the CMG actuating system. In addition, the objective func- 
tion is written to include a term insuring that CMG trajectories avoid singular 
configurations. 

The approach to momentum management studied in this chapter generalizes 
existing results in two ways: First, the spacecraft model (2.35) is used directly 
without the assumptions that accompany a Taylor series linearization, and second, 
the CMGs and singularity avoidance are treated as an integral part of the mo- 
mentum management scheme. Note that the optimal control problem formulated 
and solved in this chapter can also be regarded as an optimal maneuver problem; 
therefore, this solution procedure may be applicable to other areas of spacecraft 
control. 

The optimal momentum management problem is formulated in Section 5.1. In 
Section 5.2, a suboptimal approach is described for solving nonlinearly constrained 
optimization problems of the type posed in Section 5.1. The approach is applicable 
when the equations of motion describing the plant satisfy the external linearizabil- 
ity conditions of Chapter 3. The original optimization problem is reduced to an 
unconstrained nonlinear program in two steps. First, a linearizing tranformation is 
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applied to the system and the objective function. In this way, a new optimization 
problem is obtained in which the dynamical constraint is a linear system. Then, 
the system dynamics are approximated by a finite sum of independent functions of 
time, and a solution for the linear equations is obtained using collocation. 

In Section 5.3, the techniques of the previous section are applied to develop an 
algorithm for solving the momentum management problem stated in Section 5.1. 
This algorithm is exercised in a number of examples presented in Section 5.4. 


4.1 Formulation of the Optimal Momentum Management 
Problem 

In this section, momentum management is formulated as a quadratic optimization 
problem. First, the characteristics of a desired state trajectory are discussed for the 
model (2.35). Basically, it is desired to keep CMG momentum near some nominal 
value while keeping vehicle attitude excursions small and avoiding singular CMG 
configurations. In situations where the spacecraft initially has a significant amount 
of stored momentum, the objective of momentum management is to transfer that 
momentum out of the vehicle/CMG system. The approach to momentum manage- 
ment solved here is similar to [14], except that a nonlinear vehicle model is used. 
The optimal control problem posed at the end of the section is to minimize the 
deviation from the desired trajectory without excessive use of control effort. 

It has already been mentioned that it is undesirable to allow the CMG system to 
accumulate stored momentum along any particular axis. A large amount of CMG 
momentum projected in a particular direction implies that the ability of the CMG 
system to control the spacecraft about that axis is degraded. An exception may 
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occur if a large momentum impulse is anticipated. Under such circumstances, it 
may be desirable to transfer momentum in a particular direction to increase the 
magnitude of the impulse that can later be absorbed. These considerations can 
be used to characterize a desired trajectory for the net CMG momentum to track, 
h,Q MG {t). The desired CMG momentum may be a constant value or a ramp to a 
desired value at some point in time. Similarly, there may be a desired attitude rate 
that the vehicle should maintain or track, U7<j($). The desired net momentum can 
be expressed as 

hi(t) = h d CM0 (t) + I TEA Mt) ( 4 . 1 ) 

where all quantities are expressed in LVLH coordinates, and I TEA denotes the inertia 
matrix in LVLH of the vehicle at the torque equilibrium attitude. For applications 
such as the Space Station, the desired vehicle rate will be the orbital rate of the 
LVLH reference frame, = u> 0 . 

Certain applications may require a spacecraft to track an inertial or LVLH at- 
titude for long periods of time. Such requirements are usually specified to within a 
certain error tolerance. In other applications, characteristics of the spacecraft may 
place limits on the attainable attitudes. For example, gravity gradient torques on 
spacecraft with large asymmetries in the principal moments of inertia may make 
it difficult to maintain an attitude far from equiliburium. The requirements for 
attitude control define a desired attitude history, that the vehicle should 

maintain. 

In formulating a performance index for the CMG configuration, initial consid- 
eration was given to a measure such as the gain function [7], or a measure of rotor 
orthogonality such as (3.24). The objective function to be formulated here was 
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originally written in the latter fashion with a term of the form 

Ncmg 2 

E [(^cmo4mg)] (4-2) 

that was integrated as part of the objective function. In developing the numerical 
algorithm for application to space station momentum management, it was found 
more efficient to compute desired gimbal angle histories, based on the desired 

CMG momentum and knowledge of the optimal configuration determined by (4.2), 
rather than deal with (4.2) directly. Even though this particular term has not been 
included in the development presented here, it is not difficult to extend the following 
formulation to include it. 

The preceeding discussion illustrated that there is a natural way to define a de- 
sired state trajectory for an orbiting spacecraft equipped with momentum exhcange 
devices. Denoted a^(f), the desired trajectory can be used to define a state error 
vector 

h L (i) h%(t) 

e(t) = x(t) - x d {t) = < q L(t) - - < q%(t) ► (4.3) 

e{t) J [ e d {t) 

Assuming that there is also a penalty associated with the use of control effort, 
the cost of a particular control/state trajectory can be quantified using an integral 
quadratic performance index 

J(x,u ) = F [ e T (t) Q(t) e(t) + u T {t) R(t) u(t)] dt + e T (T) P(T) e(T) (4.4) 
where Q(t) and R(t ) denote time-dependent weighting matrices, u = 6 is the control 
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variable in (2.35), and T denotes a fixed final time. The matrix P(T) is an additional 
penalty on the error at the terminal time that can be used to drive the system to 
desired end states. The optimal momentum management problem is to minimize the 
performance index (4.4) subject to the constraint that the control/state pair satisfy 
the differential equation (2.35) with initial condition x(t = 0) = x 0 . In the next 
section, a general procedure is discussed for solving a broad class of optimization 
problems of the type posed here. 


4.2 Formulation of a Class of Optimization Problems as 
Nonlinear Programs 

In this section, an approach is developed for obtaining approximate solutions to a 
class of optimization problems which includes the optimal momentum management 
derived in the previous section. The classical approach to this type of optimization 
is to employ variational methods to obtain a set of necessary conditions that must 
be satisfied by an optimal control and its associated state and costate trajectories 
[42], [43]. The necessary conditions are in the form of a set of coupled, generally 
nonlinear, ordinary differential equations with boundary data specified at both the 
initial and final times. This two-point boundary-value problem must be solved to 
obtain an explicit solution for the optimal control. 

There are a number of difficulties that arise when one attempts to solve the 
boundary value problem for the optimal control. Because boundary data for the 
state and costate equations are specified at different times, the equations cannot 
be solved directly. In addition, different types of optimality criteria and terminal 
constraints lead to different types of boundary value problems. In general, a method 
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that works well for a particular class of equations may not be applicable to similar 
equations with different boundary data. 

A variety of gradient based methods have been developed to solve the two-point 
boundary value problems arising in optimal control theory. The ultimate conver- 
gence of these algorithms is often strongly dependent on a starting solution supplied 
by the programmer. Another difficulty is related to the fact that gradient-based 
methods can be sensitive to computational errors, and they become increasingly 
ill-behaved as the order of the problem increases. Because of these limitations, such 
techniques have not been extensively employed in autonomous control systems, and 
applications have been primarily limited to research investigations. 

Dynamic programming is another method that has been proposed for solving 
optimization problems. The difficulty with this method is that the computational 
burden increases much more rapidly than the order of the model [44]. Even though 
a number of algorithms have been proposed for resolving many of the difficult issues, 
none has proved generally reliable, and the search for methods for solving optimal 
control problems for high-order, nonlinear systems is still an active research area. 

In this section, a method is proposed for transforming a class of optimization 
problems of the form 


min 

u 

subject to 


S(x(tf) 1 tf) + / L(x,u y t)dt 

V to 

m 

X = f(x) + Y, 9i(x)ui 
»= 1 

C(x(£/), tf) = 0 
x(to) = Xq 


(4.5) 
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where S(xf,tf) is the penalty on the terminal state of the system, and L(x,u,t) is 
the penalty for a particular state/control pair. The first constraint requires that 
the state and control histories satisfy the equations of motion of the system to be 
controlled, and the constraints C(x(tj),tf) denote restrictions on the terminal state 
of the system. Limits on the magnitudes of certain state variables or controls can be 
included implicity in the objective L(x, u, t) using penalty functions. The specific 
class of problems to which the proposed method is applicable are those for which 
the vector fields /, <7i , • . . , 9m a re feedback linearizable along the optimal trajectory. 

The idea behind the development which follows is to transform (4.5) into an 
unconstrained nonlinear program. The advantage of the nonlinear programming 
formulation is that there exist robust algorithms for these problems that have guar- 
anteed convergence properties. The transformation is achieved in two steps. The 
first step uses the feedback linearization theory to transform the stated problem into 
a related problem in which the dynamical constraint is a linear system. The second 
step is to use collocation on a finite basis of functions to obtain an approximate 
solution to the linear system of equations. In applying the transformation to solve 
the optimization, care must be taken that optimal solutions remain optimal under 
the linearizing transformations, and that the collocation procedure yields a solution 
that converges to the optimal solution as the order of the approximation increases. 
These technical considerations are examined in the following two subsections. 


4.2.1 Equivalence of Optimal Control Problems 

Previously, feedback linearization has been presented as a method of transforming 
nonlinear systems into systems which exhibit linear input-output behavior. There is 
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a deeper interpretation to this problem that has not been discussed. The interpre- 
tation comes from the observation that the set of transfomations (F, a,/?) defines 
an equivalence relation on systems of differential equations of the form 

* = /(*) + £ 9i{*) v (4.6) 

X — 1 

Two systems of equations are in the same equivalence class if they are related by 
an admissible change of coordinates. Put quite simply, any system in a particular 
class can be made to behave like any other system in that class by using the appro- 
priate transformation (F, a,/?). The conditions that guarantee the existence of a 
solution to the feedback linearization problem guarantee that the equivalence class 
associated with the nonlinear system (B.52) contains one particular linear system - 
the Brunovsky canonical form (if the equivalence class contains one linear system it 
contains many). Naturally, there are different definitions of equivalence associated 
with different classes of admissible transformations [31], [32], [19], and there is also 
no reason to limit this discussion to systems which can be linearized [45]. 

Clearly, the class of transformations that have been studied for nonlinear systems 
can be used to transform (4.5) into a related problem in which the model for the 
plant is linear. Direct substitution of the functions 

X = F~\z) and u = a(F~ 1 (z)) + ^(F _1 (z)) v (4.7) 
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into the equations (4.5) yields the optimization problem 


min 

v 

subject to 


S(z{t f )) + f L(z,v,t ) 

Jt f 


dt 


z = Az + B v 
C(z(t f ),t f ) = 0 
z(to) = zo 


(4.8) 


Formally, there is no problem with this substitution; however, one would like to be 
certain that an optimal solution of (4.8) is an optimal solution of (4.5) under the 
mapping (4.7). 

The Hamiltonian function associated with the original optimization (4.5) is de- 
fined by 

H = L(x,u,t) + \ T (f(x) + G(x)u) (4.9) 

where the n additional variables A are the costate or adjoint variables. It is well 
known that an optimal solution necessarily satisfies the conditions listed below [43]. 


A 

an 

d u 

x(t 0 ) 

m) 

C(x(t/), tf ) 


f( x ) + X^,(x)u t - 

i=i 


d£ 

dx dx 



A - 



0 


X 0 

’as T dc l T 

dx dx 

L it, 


(4.10) 

(4.11) 

(4.12) 

(4.13) 

(4.14) 

(4.15) 
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The set of algebraic equations (4.12) are the stationarity conditions that determine 
the optimal control in terms of the state and costate trajectories. Equations (4.10) 
and (4.11) together with the the boundary conditions (4.13) and (4.14) form a 
two point boundary value problem that determines the optimal state and costate 
trajectories. The constant vector v denotes a set of parameters that are chosen to 
satisfy the side conditions (4.15). 

A similar set of necessary conditions exists for the problem with feedback lin- 
earized dynamics (4.8). These necessary conditions are derived from the Hamilto- 
nian 

Hi = L(z,v,t) + p T (Az + Bv) 

(4.16) 

= L[F~\z), a(F~\z)) + f3{F~\z))v, t] + p T {Az + Bv) 
and are given by 


i = Az + Bv 

(4.17) 

. T ( dL Y 

p = -A p - 1 — 1 

(4.18) 

d H 

(4.19) 

!h ~ 0 

z(t 0 ) = z 0 = F(x o) 

(4.20) 

, x \dS T dC 1 T 

(4.21) 

= aJ + " a 7, 

L t f 

= 0 

(4.22) 


In the following discussion, it is shown that if (z,v) satisfies the above conditions 
and the two systems are equivalent along z(t ), then the state-control pair (x,u) 
derived from (4.7) satisfies the necessary conditions (4.10) through (4.15). Note that 


115 



because the relations (4.7) are invertible, the converse must also hold. The necessary 
conditions are transformed explicitly in order to highlight the relationships between 
the two problems. A fairly obvious sufficient condition is also presented that only 
relies on the continuity of the linearizing transformation. This sufficient condition 
justifies solving the original optimization problem by solving the problem with linear 
state dynamics, and then transforming the solution into the original coordinates. 

The fact that the state history of the nonlinear system is related to the solution 
of the linear system should be obvious from the definition of the linearizing trans- 
formation. Direct substitution of the feedback transformation into (4.17) yields 


[D F(x)]x = [D F]{f + Ga)(x) + [D F]G0(x)v 

x = f(x) + 6r(x)[a(x) + /3(x)v] (4.23) 

= f{x) + G(x)u 


The real matrix DF denotes the derivative of the vector valued function F. Later 
results will make use of the fact that D F is nonsingular and d -~ z - = (see the 

appendix on feedback linearization for more detail). Since x(t 0 ) = F~ 1 (z(t Q )) and 
u{t) = a + /3u(J), the trajectories must be related by x = F -1 (z) on the interval 

The equation for the costate p is obtained from the derivative of the Hamiltonian 
(4.16) with respect to the state variable z. This derivative can be expanded to yield 


p T = - p t ^-{Az + Bv} - ^ 


dz 


dz 


p T —{[DF](f + Ga) + [D F]G(3v} o F~\z) 


dz 


dL 

dx 


[DF" 1 ] o F-'(z) 


dx 
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Extending the usage of the Liebnitz notation to include 

|-[D F-Ga) = ^[BF.G\a + [BF'G\~[a] 

m 

= YH DF • G V a i + [DF-G][Da] (4.26) 

1 = 1 

(4.24) can be further simplified 

f[BF] = {DF-/ + D F-Ga + DF • G/3v} 

+pT {^ [DjP ■ Ga] ~ J^ DF ' G]a } 

+? T {£lVF-GM-llVF-G]0,} - g 

= -F T l x { DF - f + T>F ^}-f X 

= -^{[sH * - i Df i s(/ + «">) - S (4 - 2?) 
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The costate equation (4.11) is recovered exactly once the substitution 


A T = p T DF 


(4.28) 


is made in the last equation. The terminal condition on A can be obtained from 
(4.21) by applying the transformation z = F~ 1 (x) and substitution of (4.28) for 
p T BF. 

The requirement (4.12) follows from (4.19) by 


dH _ dH du _ OH 
dv du dv du 


(4.29) 


Since the matrix (3 is nonsingular, multiplication on the right by its inverse yields 
the desired result, = 0. 

The preceeding argument shows that an optimal solution of the linearized prob- 
lem maps to a solution of the nonlinear problem that satisfies the required necessary 
conditions. The additional conditions which guarantee that a solution is a minimum 
- rather than a saddle point or maximum - do not transform directly to similar con- 
ditions on the equivalent problem. These conditions generally involve statements 
about second or higher-order variations of the objective function, which in turn de- 
pend on higher-order derivatives of the transformations. The only assumption that 
has been made regarding the transfomations is that both maps have nonsingular 
first derivatives. 

The difficulty with obtaining equivalent sufficient conditions is due in part to 
the fact that convexity is not preserved by the linearizing transformation. Consider, 
for example, a convex function on R that behaves like x near x = 0, and the local 
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change of coordinates x = sinz. Clearly, the composition of the convex function 
with the coordinate change cannot be convex at zero because the second derivative 
of sin z changes sign. 

Fortunately, it is only necessary to show that the transformations cannot turn 
minima into maxima or saddle-point solutions. Let J(z,v ) denote the value of 

the objective function associated with the optimal solution ( z,v ). If this solution 

is a local minimum, then there is an open neighborhood of this solution on which 
J(z,v) < J(z + 6z,v + Sv). Because the transformation (z,v) t-» (x,u) is continuous, 
there is an open neighborhood of (x,u) on which J(x,u) < J(x + 6x,u + 8u ) and the 
solution is a local minimum. Note that unless a transfomation is globally linearizing, 
uniqueness of an optimal solution is not preserved by the class of transformations 
that are being considered. 

4.2.2 Solution of Initial Value Problems by Collocation 

Functional approximation techniques have long been used in engineering and science 
to obtain solutions for boundary value problems of the form 

A(x,t) = f(x,t ) V t e D (4.30) 

B(x) = g(t ) V t G dD (4-31) 

In these equations, A denotes a differential operator on a domain D , and B is a 
constraint on the values of x imposed on the boundary dD. Such equations occur 
commonly as models for processes in heat transfer, fluid dynamics, and mechanics. 
For example, x might denote temperature, flow velocity, or strain, while t denotes a 
position in a medium. Approximate solutions to this system of equations can be ob- 
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tained by expressing the solution as a weighted sum of functions in the independent 
variable 

x*{t) = £ x { <t> { [t) (4.32) 

>=o 

The functions <f>i(t) are usually chosen to have some desired properties. Often 
they are required to satisfy the boundary conditions of the problem for all choices 
of the coefficients X{. A solution for the boundary value problem is obtained by 
substitution of (4.32) into (4.30) and choosing the X{ to minimize 

lk(*<)ll = P(*V) - /(**\0ll (4-33) 

This is a straightforward problem in functional approximation [46]. If the operators 
A and / are linear, a solution can be obtained by simply matching coefficients of 
the function inside the norm symbols. If it is required that the error residual e be 
orthogonal to the approximating functions, then the solution procedure the well- 
known Galerkin method [46], Another technique is to require that the error function 
e be identically zero on a selected set of points. This method is called collocation 
or interpolation [46]. Note that each of these techniques can be classified in some 
sense as a particular example of the Method of Weighted Residuals [47]. 

Note that the initial value problem associated with the optimization (4.5) is a 
special case of (4.30) and (4.31), and the methods described above can be applied 
to obtain a solution. These methods are not commonly applied to initial value 
problems since there are many simple forward differencing algorithms that yield 
accurate results. Also, approximation methods require many terms of the approxi- 
mating series for accurate results when the exact solution is not smooth. Experience 
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gained in this research, as well as comments by a number of authors, [48] indicate 
that when the functions involved are smooth, the approximation methods very of- 
ten give more accuracy for a given level of computational effort than differencing 
methods. 

Collocation methods are particularly simple to apply to the solution of initial 
value problems of the form [49], [50]. 

x = f(x,u) , *(<o) = x 0 (4-34) 

Approximating x(t) and u(i) using (4.32), and substitution into the previous ex- 
pression yields 

p 

X x i<M f ) ~ *p,u !,..., Up, <) (4.35) 

i=0 

p 

E*<Mo) = x ° (4.36) 

i=0 

Though it is not necessary to satisfy the the initial condition exactly, it seems 
reasonable to do so. The collocation procedure is to choose a set of times, <i, . . . ,f p , 
at which the trial functions, </>;, are independent, and then chose the coefficients 
Xi to satisfy (4.35) exactly at those times, while picking x 0 such that the initial 
condition (4.36) is satisfied. If the trial functions are complete in an appropriate 
norm, then the truncated solution converges to the exact solution of the equation 
as p increases. Using this procedure, differential equations are transformed into 
systems of algebraic equations 

p 

o = - /(®j,..., * p ,ui,..., Up, ti) 

i=0 
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(4.37) 


P 

0 = — /(®1 > • • • j ®in u l> • • • > u pi tp) 

»=0 

0 = Ys X < &(*<>) ~ *° (4.38) 

»=0 

that can be solved using readily available numerical procedures. If x is n-dimensional, 
then (4.37) with the constriant on the initial value (4.38) will consist of n(p + 1) 
equations in n(p + 1) unknowns. 

When the initial value problem is linear and x 6 R n , it is convenient to write 
the approximation (4.32) in the form 


Xq 


X(t) « * m (t)X = <f>o{t)In ••• <f> P m 


(4.39) 


Approximating the inputs in a similar fashion, u(t ) ^ $ u (t)U, the differential equa- 
tion is satified exactly at tj if 


[*.(*,-) - Afc.te)] X = B$ u (tj)u (4.40) 

If the trial functions are independent at the initial time and at the collocation 
points f p , then the requirement that the initial condition be satisfied, together 
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with the system of linear equations 


- Ai,] ((,) 


„(<i) 

• 1 

: 

X = 


[*. - A$ x \ (t p ) 


0 

B* u (t p ) \ 


(4.41) 


has a unique solution X that can be found by inverting the square matrix on the 
left. 


4.2.3 Transformation of the Optimization Problem to an Unconstrained 
Nonlinear Program 

In this section, linearizing transformations and collocation are exploited to obtain 
suboptimal solutions to the optimization problem (4.5). The basic idea is to use 
feedback linearization to transform the optimization (4.5) into (4.8), and then to 
sustitute an approximation of the form (4.32) into (4.8). The constraints become 
algebraic equations, and the integral performance measure can be evaluated using a 
suitable quadrature formula. In this way, the optimization problem is transformed 
into a nonlinear programming problem. In the absence of terminal constraints, 
the linear equations for the state approximation, Z, can be solved in terms of the 
approximation to the control V . If there are nonlinear terminal constraints, they 
can be incorporated into the objective function using Lagrange multipliers. In either 
case, the resulting programming problem is unconstrained. 

A number of variations to this basic approach have already appeared in the lit- 
erature. The methods differ primarily in the way in which the original problem is 
approximated, and in the numerical method used to solve the problem. Solutions 
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have been approximated using orthogonal polynomials, splines, and Fourier series. 
In some cases, both the control and the state variables are parameterized, and a 
solution is obtained by solving the resulting constrained nonlinear program. In [51], 
the powerful method of successive quadratic programming [52] is applied directly 
to the constrained nonlinear program. Because the authors use cubic splines to 
approximate the time functions on intervals, it is necessary to satisfy rather cum- 
bersome compatibility conditions at the end points. These conditions require that 
four equations be solved for each collocation point. Another procedure incorporates 
the constraint into the objective function using Lagrange multipliers [48]. The mul- 
tipliers effectively double the size of the programming problem, which may already 
be quite large (~ np variables requiring ~ ( np ) 2 storage). 

Other techniques parameterize only the state variables or only the control vari- 
ables and solve the resulting algebraic equations for the control or state history, 
respectively. In [27], truncated Fourier series are used to parameterize the state 
trajectory for a class of systems of the form 

F(x,x,x) = u (4.42) 

The simple form of the system equations makes it possible to eliminate the control 
from the optimization problem. Two difficulties with this approach are that com- 
pleteness of the approximating series becomes a delicate issue when parameterizing 
the state history, and a number of technical difficulties arise when the equations are 
over- or under-determined in the control variable (especially under-determined). 

Note that the choice of approximating polynomials is not critical to the solution 
of the initial value problem by the collocation method. Any set of functions which 
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are independent and complete in an appropriate sense will work. In the published 
literature, Chebyahev polynomials are used quite extensively [49], [48] 1 (see [53] for a 
detailed discussion of the approach using these functions). In certain circumstances, 
the particular problem at hand may suggest a choice of functions that should yield 
accurate results [54] (boundary value problems), or perhaps simplifies the compu- 
tations. The polynomials that have been chosen for use here are the interpolating 
polynomials given by 

£i(t) = n l _ ! r z = ( 4 - 43 ) 

i=0 li l 3 

i ^ j 

and the polynomial approximation to the state is given by 

x{t) « Y, x k £ k (t) (4.44) 

k=0 

It is assumed that < 0 is the initial time and the coefficient x 0 is chosen to satisfy 
the constraint at the initial time, x 0 = z(f 0 )- The x± for i > 1 are to be chosen by 
collocation at the times for l <i < p, which are chosen to be the zeros of a p th - 
order Legendre polynomial defined on the interval [<o,i/]. This choice of collocation 
times was made because of the fact that the product defining the interpolating 
polynomials ranged from 1 to p, and these zeros were used for the times tj y then the 
resulting set of functions would be orthogonal on the interval [to,tf]. In addition, 
because the zeros are more closely spaced at the end points of the interval, the 
approximation is better behaved near the end points of the interval (see [46] and 
the discussion of the Runge phenomenon). 

l Note that the dates on these references are 1956 and 1988, respectively. 
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Note that it is also possible to include the final time, tf , as one of the interpola- 
tion times. In some cases, this was found to improve the appearance of the plotted 
solutions by making the last segment of the time histories better behaved. The 
effect on the value of the objective function was negligible, provided the choice of p 
was sufficiently large. Including the final time as a collocation point was found to 
be undesirable when terminal costs were included in the cost function. This tended 
to make the gradient of the objective function very steep in a particular direction, 
and adversely affected convergence of the algorithm. 

The advantage of the interpolating polynomials derives from the fact that £i(tj) = 
0 unless i = j, which implies that the value of the approximation at a collocation 
point ti depends only on x;. The collocated approximations to state dependent 
functions, such as the cost functional and the linearizing transformations, can be 
evaluated by considering only one finite set of data, ( z;,u; ), at a time. The lineariz- 
ing transformation as well as the approximate solution of the nonlinear dynamics 
can be evaluated using 


X; 


F-\zi) 


«(*(<)) 

(3{x(t)) 

x(t) 

u(t) 


X «(*i)4(0 

t=0 

X 0(*O 4(0 

i = 0 

E f'teW) 

i = 0 

X ( a (*0 + 0(*iM) 4(0 


i=0 


( 4 . 45 ) 

( 4 . 46 ) 

( 4 . 47 ) 

( 4 . 48 ) 

( 4 . 49 ) 

( 4 . 50 ) 
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A similar set of expressions exist to invert these transformations. Computing the 
approximation of the transformation (z,v) (x,n) reduces to solving p indepen- 

dent sets of n equations in n unknowns, as opposed to one system of up equations 
that would result if another set of polynomials had been chosen 2 . 

It is now simple to state the nonlinear program which yields the approximate 
solution to (4.5). Letting L{ = L(xi,Ui } ti), the integral performance index is ap- 
proximated by 

L(x, u, t)dt ^ y; Lj 

= 5Z (jf ti (t)dtj 

= ii (4.5i) 

i=0 

(4.52) 

The integral of £i(t) has been evaluated numerically for p < 20 and is positive for all 
i. A formal proof that this is true for all p is not possible since a closed-form solution 
for the interpolation times is not available; however, inspection of the functions 
involved suggests that this is true for larger p . The terminal cost can be evaluated 
as S{z{tf) } tf), and the terminal time must satisfy the constraints C(z(tj),tf). Since 
the Z{ are derived from a linear matrix equation of the form Z = Zo + D V y the 
approximate cost and the terminal constraints can be expressed entirely in terms 
of the mp - dimensional vector V. Letting A denote a vector of Lagrange multipliers, 
the nonlinear programming problem to be solved for a suboptimal solution to (4.5) 

2 The initial time is not counted with the others 
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IS 


Y, L ( z i,Vi,ti)£i + S(z(t f ),t f ) + X T C(z(t f ), t t ) (4.53) 

A, V i = o 

In the following section, the optimal momentum management problem is expressed 
in this form, and a computational algorithm is outlined to solve the nonlinear pro- 
gram. 


4.3 Solution of the Optimal Momentum Management Prob- 
lem 


In this section, the technique of the previous section is applied to obtain subopti- 
mal solutions to the optimal momentum management problem (4.4). The integral 
quadratic performance measure can be converted to an algebraic function using the 
polynomial approximation. Letting e(t) « ]T?=o e»^(<) = £(t)-E\ and considering 
only one of the terms in the performance index (4.4), the criteria can be written 
as 



e(t) T Q e(t) dt 
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■ e P (t)i n ]dt\ 


(4.55) 
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(4.56) 

(4.57) 


= E t QQE 
= **(Z)QQV e (Z) 


Because E is a vector of constants, it can be taken outside the integral in (4.54). 
If Q is a function of i, then this equation can be integrated during the initialization 
of the optimization algorithm to obtain a constant matrix. With Q assumed to be 
constant, it can also be extracted from under the integral by rewriting the matrix 
product. In (4.55), the integral depends only on products of the approximating 
polynomials, and a closed-form solution is easily obtained for each fixed value of p. 
The fact that only one of the interpolating polynomials is non-zero at each of the 
collocation points can be used to express ^S e {Z) in the form 


X - X d = * e (Z) 


F-\z o) 

< : 


[ F~ l {z p ) 


(4.58) 


An expression similar to (4.57) can be derived for the integral of the term quadratic 
in the control in (4.4). This expression has the form 



u T (t)Ru{t)dt ^ *Z(V)RRV U (V) 


(4.59) 


where 


*»(v) = 


a(F- 1 (z 0 )) + PiF-'izo)) v 0 


a(F-\z p )) + 0{F~\z p )) Vp+1 


(4.60) 


129 



Finally, the cost for the terminal error is simply expressed as 


e r (T) P e(T) = E T C t (T) P> P> E C(T) (4.61) 

Note that each of the terms in the cost functional has been reduced to a nonlinear 
quadratic form which can be expressed as 

n(p+2)+m(p+l) 

min F(V) = min fi'fi (4.62) 

V V i=i 

where the / t - are the components of a vector valued function /. 

Figure 4.1 shows a flow chart of the computational algorithm that has been 
implemented to solve the optimal momentum management problem. The initializa- 
tion routine reads data from an external file and computes an initial guess for the 
optimal solution. Along with initial state data, a data file is read that includes the 
integrated polynomial products £i(t) ij(t ), and the matrix W that solves the collo- 
cation problem for the state history Z = W(Z 0 + V). The initial guess is obtained 
by initializing the state trajectory of the spacecraft so that the vehicle follows a 
torque equilibrium attitude, the momentum remains fixed in inertial coordinates, 
and the CMG gimbals rotate to hold the rotors fixed in inertial space. This solu- 
tion is transformed to an equivalent state history for the linearized system, and the 
control producing the best approximation to this solution in a least-squares sense 
is the starting solution. 

The main sequencing routine then calls a subroutine DUNLSJ 3 that solves the 
nonlinear quadratic programming problem using the Lavenberg-Marquardt algo- 

3 IMSL subroutine library, version 10.0 
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Figure 4.1: Flow Chart of Optimal Momentum Mangagement Algorithm. 
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rithm. This routine calls two subroutines that compute the components of the 
function / in (4.62) as well as the Jacobian of /. The function is evaluated in 
several steps. First, a solution to the linear system dynamics is obtained, and this 
solution is transformed to obtain the solution in terms of the approximation to the 
nonlinear system dynamics X. Then the feedback transformation is computed, and 
the vector U is evaluated. The components of the function f are evaluated and exe- 
cution is returned to the optimization subroutine. As mentioned above, the routine 
DUNLSJ also calls a subroutine that evaluates the Jacobian of the function /. The 
solution is assumed to have converged when one of the following three conditions 
holds: i) The value of the function F(V) is below a certain threshold, ii) the gradient 
of the function is below a certain threshold, and iii) three successive iterations have 
failed to improve the objective function by more than one percent. In the following 
sections, results are presented that have been generated using this program. 
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4.4 Sample Results of Optimal Momentum Management 


In this section, a number of examples are presented to illustrate the benefits obtained 
by solving the optimal momentum management problem. The primary advantage 
of the proposed approach is that issues relating to attitude control and momentum 
management are treated in a combined framework which includes a nonlinear vehicle 
model. 

A number of continuous time, linear optimal momentum management schemes 
have appeared in the literature [13], [14], [55]. These schemes assume that both 
attitude and momentum excursions will remain small. Assuming that the spacecraft 
is to track a torque equilibrium attitude, the former assumption is probably valid. 
The latter assumption is generally not justified; results presented in this section as 
well as [13], [14] themselves indicate variations in net vehicle momentum of the same 
order of magnitude as the nominal values. 

Another approach to momentum management involves computing a specified 
amount of momentum that should be transferred during an orbit. The desired 
momentum transfer is achieved by manipulating the attitude dependent gravity 
gradient torque [15], [55]. The momentum to be transferred is computed by adding 
the amount of momentum currently stored in the CMGs to the momentum that is 
expected to accumulate over the next orbit. These methods consider neither the 
constraints imposed by vehicle dynamics nor advantageous ways to exploit vehicle 
dynamics. Both cited references neglect the fact that the inertial orientation of 
the gravity gradient torque is changing continuously. Since these methods are also 
based on a linearized model, they are not applicable when there are large attitude 
excursions. 
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One of the most important features of the optimal momentum management 
algorithm is the way in which the solutions exploit the total vehicle momentum 
iself for regulating both momentum and attitude. This feature is highlighted in 
the first set of examples in which optimal solutions are obtained with the space- 
craft subjected to various environmental torques. In the first example, only the 
gravity gradient torque is considered. The next example includes a constant bias 
torque about the yaw axis, while in the third example, the torque was chosen to be 
representative of aerodynamic disturbances the space station is likely to encounter 
on-orbit. The final two examples show solutions in which the CMGs are initialized 
with net momentum storage in the pitch and yaw directions. 

The spacecraft model used in the examples presented in this section is the same 
as that employed in the steering examples of the previous chapter. The vehicle 
inertia properties are listed in Table 3.2 and represent data for the rephased dual- 
keel space station, after assembly sequence fight 7 [40]. This vehicle is illustrated 
schematically in Figure 3.2, and the CMG mounting configuration is shown in Fig- 
ure 3.3. 

In each of the examples to be discussed, the vehicle is in a circular orbit at an 
altitude of 300 fcm, and the orbital period is 90 minutes . The momentum manage- 
ment problem is solved for one orbital period. The desired solution is computed by 
assuming that the vehicle is to track the torque equilibrium attitude, and that the 
CMGs are to track a zero momentum state. For the vehicle data used in this ex- 
ample, the gravity gradient torque required the equilibrium attitude to be pitched 
almost -17 degrees from the local vertical. The desired yaw attitude is assumed 
to be zero. In order to track an orientation that is fixed in the LVLH reference 
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frame, the vehicle must rotate with respect to an inertial frame at a rate equal 
to the LVLH or orbital rate. As a consequence, the momentum of the spacecraft 
will not be zero, and because of assymetries in the mass distribution, the vehicle’s 
angular momentum need not align with the LVLH rate. The desired momentum at 
the initial time is computed by determining the momentum of the vehicle as if it 
were fixed at the desired attitude, and rotating with the LVLH frame. The desired 
momentum trajectory is then computed such that the momentum at the initial time 
remains constant in inertial coordinates as the LVLH frame rotates. The desired 
gimbal angle history is computed by assuming that the CMGs are to track the zero 
momentum configuration initialized with gimbal angles (in radians) 
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( 4 . 63 ) 


The outer gimbal angles are assumed to rotate at the orbital rate. Because the 
CMGs are mounted with outer gimbal axes along the negative vehicle pitch axis 
(which aligns with the orbital rate), the desired trajectory causes the individual 
rotors to remain fixed in inertial coordinates. 

It should be remarked that a number of the optimal solutions exhibit trends 
near the terminal time that seem to be inconsistent with optimal behavior. In 
particular, there can be large excursions in the state errors and controls near the 
final time. These appear to be due to properties of the polynomials being employed, 
and are probably related to the Runge phenomenon that occurs in polynomial 
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approximation [46]. For a given optimization problem, the precise behavior at the 
final time was seen to change noticably as the order of the approximation increased. 
This occured even though the magnitude of the objective function appeared to have 
converged quite close to the optimal value, and differences between solutions were 
hardly discernable, except near the endpoint. Including an interpolation point at 
the final time improved the behavior there, but slowed convergence of the objective 
function. There are two potential remedies for this difficulty; the first is to include 
a penalty on the final error, and the second is to extend the optimization beyond 
the orbital period to obtain useful solutions for the time interval of interest. 


4.4.1 Baseline Solutions Illustrating Required Computational Time 

Before considering specific examples, a generic set of runs are investigated to illus- 
trate how the amount of processor time required to obtain solutions to the nonlinear 
program varies as the order of the approximation, p , was increased. In these cases, 
it is assumed that only gravity gradient torques act on the vehicle, and that the 
initial attitude, as well as the desired attitude, correspond to the equilibrium orien- 
tation with zero yaw. The desired momentum trajectory is computed by requiring 
that the momentum should remain constant in inertial coordinates. The desired 
gimbal trajectory was generated by rotating the CMGs about their outer gimbal 
axes at the orbital rate. The penalty on the terminal error, P, is set to zero, the 
weighting on the control, J2, is a diagonal matrix of 10 +4 , and the error weighting, 
Q , is a diagonal matrix of 10 -4 for the momentum errors, 10 +4 for the attitude 
errors, and 10 +2 for the gimbal angle errors. These weightings were chosen so that 
the contribution to the total cost associated with the attitude error, momentum 
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error, gimbal angle error, and control effort are of approximately the same order of 
magnitude in the absence of disturbance torques. 

Table 4.1 shows the variation of the required computation time and the value 
of the objective function as the order of the approximation is increased. The CPU 
times are for a Digital Equipment Corporation VAX 8600 4 . Note that the solutions 


p 

CPU time 

Objective 

% Improvement 

3 

9.72 

3174 

- 

4 

17.94 

2781 

12.5 

5 

21.74 

1603 

42.4 

6 

27.79 

1603 

0.0 

7 

90.71 

1588 

0.9 

8 

80.61 

1586 

0.1 

9 

124.63 

1584 

0.1 


Table 4.1: Nominal CPU times require to solve optimal momentum management 
problem. 


for p = 5,6 have the same value of the objective function, even though the value 
should decrease with increasing p. This is probably due to the fact that for p < 6, 
the algorithm terminates when the gradient of the objective crosses a threshold. 
For larger values, it terminates when three iterations fail to reduce the objective 
by more than one percent. The CPU time required for the case with p = 7 to 
converge is clearly inconsistent with what should have been expected; although this 
was observed in another situation, no satisfactory explanation has been found. 

The approximation using p = 8 has been used most extensively, and required 
computation times were typically less than the 80 sec. indicated in Table 4.1. It 
was generally more common to see times in the range 40 - 60 sec. for p = 8, and 

Conversion factor to equivalent VAX 11/785 CPU seconds is ^ 3.5 
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the numbers shown in the table might be regarded as approximate upper bounds. 
It is interesting to note that the required convergence time was not very dependent 
on how near the initial guess was to the optimal solution. The number of iterations 
tended to increase with p, but was fairly consistent for a fixed p, and almost never 
exceeded 10. Because of a number of large matrix computations that must be 
performed, the percent increase in CPU time is larger than the percent increase in 
P* 

Extensive experience using the computer program developed to solve the optimal 
momentum management program has shown it to be very robust. In contrast 
to methods that solve the two-point boundary value problems arising in otpimal 
control, the algorithm developed here was found to be quite insesitive to the choice 
of starting conditions. In all the results presented here, the same initial guess is 
used to start the omptimization. When the algorithm did fail to converge, it was 
generally necessary to vary only a parameter that limited the amount the solution 
could be changed for each iteration. While it was desired to keep this limit as large 
as possible, situations occured where the algorithm oscillated about the optimum 
or an initial step was so large that the linearizing transformations could not be 
computed. 
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4.4.2 Momentum Management with Gravity Gradient Torque Only 


In this example, the momentum management problem is solved assuming that grav- 
ity gradient torques are the only external torques acting on the vehicle. The results 
in this example show how the optimal solutions use CMG momentum transfer and 
gravity gradient torques to offset the Euler torque acting on the vehicle. The ex- 
ample also provides a baseline for comparison with later results. 

Because state excursions are small, it is reasonable to assume that a linearized 
model would yield approximately the same results shown here. Since the CMGs 
are initialized in a zero momentum state, and there are no external torques, a 
momentum manager like [15] based on a momentum dumping scheme would choose 
to follow the gravity gradient equilibrium attitude. Since the attitude control law 
of a vehicle using this type of momentum manager should compensate for Euler 
coupling, these results will also be similar. Note that both of these schemes assume 
that the steering law is causing the outer gimbals to rotate at the orbital rate. 
While a well designed steering law will have this feature, not all steering laws can 
guarantee this behavior, and momentum may accumulate in the CMG system. 

For this example, the vehicle is initialized at the torque equilibrium attitude with 
the CMGs in the configuration (4.63). The total LVLH momentum components are 
plotted in Figure 4.2. This figure shows that the momentum rotates in the LVLH 
frame so that it remains fixed in inertial coordinates. Because the orbital rate aligns 
with the pitch component of momentum, this component remains constant. Since 
the initial momentum had a small component along the negative roll direction in 
the LVLH frame, the roll and yaw components appear as out-of-phase sine-waves 
due to the rotation of the LVLH reference frame along the orbit. 
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A plot of CMG momentum in LVLH coordinates is shown in Figure 4.3. Initially, 
the Euler torque due to rotation of the LVLH frame (u; x Iu>) is on the order of -0-2 
ft-lb in the LVLH yaw direction, with a very small component in the positive roll 
direction. The Euler torque is caused by the rotation of the spacecraft’s angular 
momentum in inertial coordinates. The spacecraft momentum is initially 287 ft-lb 
in the negative roll direction in LVLH. To remain fixed in inertial coordinates, the 
momentum must rotate with respect to LVLH at a rate equal to the orbit rate. Since 
Euler coupling is not considered in computing the initial or desired solutions, the 
CMGs must dynamically compensate for its influence during the orbit. Figure 4.3 
indicates that the CMGs transfer momentum at a rate approximately equal to the 
rate of change of the vector Ioj in LVLH coordinates. This produces a torque on 
the spacecraft approximately opposite to the Euler torque. 

The plot of attitude error, Figure 4.4, indicates that the optimal solution favors 
a positive displacement about yaw and roll. The displacement about roll causes a 
small positive bias in the gravity gradient torque about that axis. Because the de- 
sired attitude is gravity gradient unstable in both pitch and roll, positive deviations 
from the equilibrium result in positive gravity gradient torques. The effect of the 
gravity gradient torque together with the momentum associated with the space- 
craft rotation causes the errors in total momentum seen in Figure 4.5. Because 
the gravity gradient torque is nearly constant in LVLH coordinates, the rotation of 
the LVLH frame yields variations in the yaw and roll components of total momen- 
tum that appear approximately as out-of-phase sine-waves. The momentum errors 
that appear in this figure cause the net momentum to become more aligned with 
the angular rate vector. This produces a reduction in the Euler coupling and a 
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corresponding reduction in CMG control effort. 

Note that some of the time histories in the figures appear to diverge near as 
the solutions approach the final time. This behavior is due to the nature of the 
approximating functions, and was discussed in some detail in the introduction to 
the examples. The behavior is also evident in the other examples presented in this 
chapter. 

When a constant bias is applied in the roll or pitch directions, the optimal 
solutions are virtually identical to the solutions shown in this example. The reason 
for the similarity is that the optimal momentum manager adjusts the attitude about 
these axes so that the gravity gradient torque opposes the external torque on the 
vehicle. Since the gravity gradient torque is orthogonal to the yaw direction in 
LVLH, constant torques in this direction yield solutions that are noticably different 
as seen in the next example. 
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MOMENTUM ( FT-LB-SEC ) MOMENTUM ( FT-LB-5EC ) 

-400.0200.0 0.0 200.0 400.0 600.0 -1000.0.0 1000.02000.03000.04000.05000.0 


TOTAL MOMENTUM. IN LVLH 



0.0 500.0 1000.0 1500.0 2000.0 2500.0 3000.0 3500.0 4000.0 4500.0 5000.0 5500.0 


TIME (SECONDS) 


Figure 4.2: LVLH Momentum: Gravity Gradient Torque Only (513). 

TOTAL CMG MOMENTUM IN LVLH 



0.0 500.0 1000.0 1500.0 2000.0 2500.0 3000.0 3500.0 4000.0 4500.0 5000.0 5500.0 

TIME (SECONDS) 


Figure 4.3: CMG Momentum in LVLH Frame: Gravity Gradient Torque Only (513). 
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MOMENTUM ( FT-LB-SEC ) ATT I TUBE *10 

- 15 . 0 - 10.0 - 5.0 0.0 5.0 10.0 - 5.0 0.0 5.0 10.0 15.0 20.0 25.0 


RTTITUDE ERROR QUATERNION urt LVLH 



TIME (SECONDS) 


Figure 4.4: Attitude Errors: Gravity Gradient Torque Only (513). 


MOMENTUM ERROR IN LVLH 



TIME (SECONDS! 


Figure 4.5: LVLH Momentum Error: Gravity Gradient Torque Only (513). 
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4.4.3 Momentum Management in the Presence of a Yaw Bias Torque 


In this example, a constant torque of 0.1 ft- lb is applied along the yaw direction in 
the LVLH reference frame. Torques are applied in LVLH coordinates rather that 
body coordinates because the momentum is expressed in this reference frame, and 
the results are more easily interpretted when it is not necessary to consider the 
influence of a nonzero attitude. The initial data and error penalties are the same 
as in the previous example, and the desired data are computed as described in the 
introduction to the examples. 

Figure 4.6 shows that a momentum error of about 100 ft-lb-sec quickly accu- 
mulates in the positive roll direction, along with a smaller error in the negative 
yaw direction. The reasons for this momentum transfer can be understood by ex- 
amining the rate of change of momentum in inertial coordinates. If the effect of 
the external torque were not compensated for, the torque would cause the net mo- 
mentum of the vehicle system to change at a rate of 0.1 ft-lb in the positive yaw 
direction of the LVLH frame. The momentum transfer in the roll direction couples 
with the orbital rate to produce a momentum change in the negative yaw direc- 
tion that approximately cancels the momentum transfer due to the external torque 
(iv 0 x h L « —0.001 • 100 « —0.1). The necessary momentum transfer is produced 
by a positive deviation from the torque equilibrium attitude that can be seen in 
Figure 4.7. 

Figure 4.7 also shows an attitude error accumulateing in the yaw direction, but 
is nulled before the final time. Since the error eventually stops increasing, the 
CMGs must be applying torques to the spacecraft to limit the attitude excursion. 
The influence of the CMGs can be understood from Figure 4.8 which shows the 
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momentum transfer in body coordinates. Initially, The CMGs compensate for the 
external torque by transferring momentum into the negative roll and yaw directions. 
The momentum transferred in the negative roll direction produces the positive 
attitude error in yaw. Because the vehicle is unstable at the desired attitude, this 
error results in a positive torque in the LVLH frame that produces the momentum 
error seen in Figure 4.6. The momentum transfer in yaw produces a small positive 
displacement in yaw. 

The attitude rate changes that occur at the beginning of the solution result in 
changes in the vehicle momentum approximately equal to the momentum transfer 
from the CMGs. Note that even a change of 500 ft-lb-sec is small compared to the 
total momentum of the orbiting vehicle (about 4000 ft-lb-sec), the net change in 
the pitch-yaw plane is about twice the nominal value. It is the momentum that lies 
in the plane perpendicular to the pitch axis that determines the amount of Euler 
coupling. A momentum management scheme based on a tangent linearization of 
the dynamical model would underestimate this torque by a factor of two. 
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ATTITUDE *10* MOMENTUM ( FT-LB-SEC ) 

0 1.0 2.0 3.0 4.0 — -100.0 -50.0 0.0 50.0 100.0 
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TIME (SECONDS) 


igure 4.6: LVLH Momentum Error: 0.1 ft-lb Torque about Yaw (566). 


RTTITUDE ERROR QURTERNION wrt LVLH 
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Figure 4.7: Attitude Errors: 0.1 ft-lb Torque about Yaw (566). 
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TOTRL CMG MOMENTUM IN LVLH 



Figure 4.8: CMG Momentum in Body Frame: 0.1 ft-lb Torque about Yaw (566). 
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4.4.4 Momentum Management in the Presence of an Aerodynamic Torque 


In this section, the torque applied to the vehicle is chosen to be representative of 
aerodynamic torques that the vehicle is likely to encounter on-orbit. The applied 
torque is shown in Figure 4.9 and was chosen to be similar to a plot in [15] which 
was generated by a computational model of the fully erected space station. The 
sinusoidal variation in the torque is due to the variation in atmospheric density 
as the spacecraft orbits on the night and then the dayside of the Earth. More 
complicated behavior can be expected if effects such as solar panel rotation are 
taken into account. 

The attitude error is shown in Figure 4.10. As the external torque about the 
pitch axis becomes less negative, the pitch component of the equilibrium attitude 
becomes more negative. The negative attitude error about this axis produces a neg- 
ative torque on the vehicle that aids in tracking the equilibrium attitude and reduces 
momentum requirements on the CMGs. The attitude error about roll exploits grav- 
ity torques to cause the transfer of momentum in the negative roll direction that can 
be seen in Figure 4.11. In the previous example it was explained how momentum 
transfer in the pitch direction aids in compensating for torque disturbances in the 
yaw direction. The CMG momentum transfer in LVLH is shown in Figure 4.12. The 
momentum transfer at the starting time in the pitch and roll directions produces 
negative torques on the vehicle and yields the resulting attitude changes. Note once 
again that the momentum excursions are large enough that a linearized model will 
miscalculate the Euler coupling by a factor of between four and eight. 
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EXTERNAL TORQUE 
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Figure 4.11: Momentum Errors: Aerodynamic Torque (606). 
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Figure 4.12: CMG Momentum Transfer: Aerodynamic Torque (606). 
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4.4.5 Momentum Management with Initial Momentum Error in Pitch 

In this example, the CMGs are initialized with a net momentum storage along the 
pitch axis (i.e.- parallel to the outer gimbal axes). Except for gravity, there are no 
other external torques acting on the vehicle. The initial gimbal angles are all - 0.1 
radians , and all other data are the same as in the previous examples. 

Figure 4.13 shows that the initial momentum error is 2000 ft-lb in the pitch 
direction, and that the momentum is smoothly removed by the end of the orbit. The 
momentum error reduction is achieved by manipulating the gravity gradient torque. 
The large pitch attitude error apparent in Figure 4.14 produces a negative gravity 
torque on the vehicle. Figure 4.15 shows that the CMGs provide the torques required 
to produce this attitude error, and then to prevent the attitude from diverging due 
to the unstable nature of the equilibrium. 

Inspection of Figure 4.15 indicates that the CMGs exchange about 5000 ft-lb- 
sec of momentum with the spacecraft along the pitch axis. As a consequence, the 
momentum stored in the rotating spacecraft, as well as the Euler torque on the 
spacecraft, more than double. Once again, a linearized model could not account for 
these effects without imposing restrictions on the allowable vehicle rates. 
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MOMENTUM ERROR IN LVLH 
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Figure 4.13: Momentum Errors: Initial Momentum Error in Pitch (569). 


RTTITUDE ERROR QURTERNION wrt LVLH 



TIME (SECONDS) 


Figure 4.14: Attitude Errors: Initial Momentum Error in Pitch (569). 


152 



TOTAL CMG MOMENTUM IN LVLH 



TIME (SECONDS) 


Figure 4.15: CMG Momentum Transfer: Initial Momentum Error in Pitch (569). 
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4.4,0 Momentum Management with Initial Momentum Error in Yaw 

The problem solved in this example is the same as the previous example, except that 
the CMGs are initialized with net momentum storage in the yaw direction of the 
LVLH reference frame. The momentum error in yaw is obtained by initializing the 
outer gimbal angles of CMG#2 and CMG#3 at ^ — 0.1 and ^+ 0 - 1 ? respectively. All 
other data are initialized as in the first example where there were neither external 
torques nor initial momentum errors. 

Figure 4.16 shows that the optimal solution adequately removes the momentum 
error that was present at the initial time. Note that there is initially an increase in 
the momentum error in the pitch direction. The reason for this can be understood 
from equation (2.56), which shows that the rate of change of momentum in the 
yaw direction is equal to the product of the orbital rate with the momentum in 
the pitch direction. The attitude errors and the CMG momentum transfer shown 
in Figures 4.17 and 4.18, respectively, show that there is a significant amount of 
momentum exchange occurring, even though the rate-of-change of net momentum 
is quite smooth. Note also the correlation between the attitude error about roll, 
and the rate of net momentum change, which is due to the gravity torque. 
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Figure 4.16: Momentum Errors: Initial Momentum Error in Yaw (570). 


RTTITUDE ERROR QURTERNION wrt LVLH 



TIME (SECONDS) 


Figure 4.17: Attitude Errors: Initial Momentum Error in Yaw (570). 
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Figure 4.18: CMG Momentum Transfer: Initial Momentum Error in Yaw (570). 
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4.5 A Conceptual Spacecraft Control System Based on Feed- 
back Linearization 

In this section, an attitude control system is described that incorporates the opti- 
mal momentum management scheme and the CMG steering algorithm described in 
Chapter 4 into a control system that regulates not only the vehicle attitude, but 
also the stored momentum and the CMG configuration. The schematic structure of 
the system is shown in Figure 4.19. Functionally, the system can be seperated into 
a command generator, a feedback control system, and a plant whose dynamics are 
to be controlled. 

In the case of interest here, the plant is a spacecraft equipped with CMGs. 
The command generator consists of the optimal momentum management scheme 
described in the previous chapter. The optimal momentum management problem 
is solved to obtain a state history for the vehicle to track over an orbit. Closed- 
loop control is required to regulate the actual spacecraft state to remain near the 
commanded state in the presence of uncertain environmental torques. The design 
of these feedback control laws has not been considered in this thesis. If errors in the 
environmental disturbances become so large that the control system has difficulty 
tracking the desired solution, the momentum manager can be called to generate a 
new solution that also desaturates any excess momentum stored in the system. 

The functions of attitude control and momentum control have been shown sep- 
arately in Figure 4.19. One reason for this is that it seems reasonable to consider 
momentum control as being subordinate to attitude control in typical on-orbit appli- 
cations. Also, the control authority over attitude dynamics provided by the CMGs 
is likely to be greater than the authority gravity gradient torques have over momen- 
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turn, implying that the two control systems may have widely differing bandwidths. 
Furthermore, there are a number of other actuation devices that may be used to 
control momentum that have not been considered here. An evaluation of the precise 
function of the feedback control systems along with procedures for obtaining useful 
designs has been left for future investigation. 

Note that the approximation based steering algorithm has also been included in 
the control system shown in Figure 4.19. In this case, the approximation procedure 
outlined in Chapter 3 should be formulated using the affine control parameterization 
for the system model that includes the external torque, and the desired rates are 
computed by the optimal momentum manager. The approximation algorithm is 
still needed because there will always be distrubances that are not accounted for in 
the optimal solution. Even though the approximation algorithm should continue to 
monitor the CMGs for undesirable configurations, its primary function here is to 
distribute disturbances among the available CMGs so that the desired trajectory 
can be tracked. 
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Figure 4.19: Schematic of an Attitude/Momentum Control System for a Spacecraft 
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Chapter 5 

Summary and Recommendations for Future 
Research 

The goal of this thesis was to develop an appoach to the control of spacecraft with 
CMGs in which the vehicle and the CMGs are treated as a combined system. In 
general, the equations for such a system are nonlinear and highly coupled. Feedback 
linearization techniques were used to transform appropriate nonlinear spacecraft 
models into equivalent linear systems. These linearized models were then used to 
develop practical schemes for solving the CMG steering problem and the optimal 
momentum management problem. 

In Chapter 2, two spacecraft models were derived along with the state and input 
transformations that linearize them. The difference between these two models was 
that one assumed the spacecraft was an isolated system with inertially constant 
net momentum. This model included state equations for the vehicle attitude and 
the CMG gimbal angles. The second model incorporated the effect of the attitude 
dependent gravity gradient torque, and provided for the accomodation of arbitrary 
external torques. In addition to equations for the vehicle attitude and CMG state, 
this model included the net momentum as a state variable. The chapter concluded 
with the development of transformations that linearized both models. 

In Chapter 3, an approach to CMG steering was developed utilizing the trans- 
formation that linearized the dynamics of the first spacecraft model. Although the 
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system to be controlled had been made linear, the transformations themselves did 
not address certain critical issues relating to the steering problem. For systems in 
which there are redundant actuators, it was pointed out that the linearizing trans- 
formations effectively eliminate the redundancy. This is particularly undesireable 
when redundancy has been included to increase the flexibility of the system to be 
controlled. To remedy the situation, a method was developed to distribute control 
effort among redundant actuators in a feedback linearized system, and then applied 
to the spacecraft with CMGs. 

The rtormed approximation algorithm developed in Chapter 3 for CMG steering 
was shown to have a number of advantages over existing algorithms. With this 
method, the problem of CMG steering law design involves defining how the CMGs 
should behave in a given situation, and then chosing a set of gimbal rates that pro- 
duces this behavior. Conceptually, this problem is simpler than the conventional 
method of determining a function whose gradient defines the desired behavior. Us- 
ing this approach, it is possible to develop a steering law that can recognize problem 
situations, and take specific action to remedy the difficuly. Another important fea- 
ture of the approach is that it incorporates explicit bounds on the rates that can 
be produced by the gimbal drives (reflecting actual hardware limitations). With 
one exception, existing approaches scale the commanded gimbal rates to obtain a 
set of rates that do not violate physical constraints, thus failing to produce the re- 
quired torque on the vehicle. A simple modification to the proposed algorithm was 
described that guaranteed the steering law would not fail when forced to operate 
near a singularity. 

Examples were presented to demonstrate the behavior of the proposed steering 
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algorithm performing a variety of attitude control functions. The first three ex- 
amples illustrated the characteristics of the gimbal angle trajectories generated by 
the three normed approximation algorithms. The 1-norm and oo-norm algorithms 
were found to exhibit undesirable oscillations in the present application. Because 
the 2- norm algorithm generally yielded smooth, well behaved solutions, it was used 
in the remaining examples. The fourth example demonstrated a situation in which 
the steering algorithm was able to deal explicitly with a problem situation that con- 
ventional gradient based methods would fail to recognize. The next two examples 
showed the algorithm redistributing the CMG rotors when they are initialized in 
poor configurations. In one of these examples, the algorithm was initialized in a 
singular configuration. The final example of this chapter showed the response of 
the vehicle/CMG system when a large onboard disturbance occurs. The distur- 
bance was large enough to saturate several of the CMGs in rate, and the example 
demonstrates the advantage of enforcing rate limits on the individual gimbals. 

In Chapter 4, an approach was developed to compute suboptimal trajectories for 
an optimal momentum management problem. The method used feedback linearizing 
transformations to restate an optimal control problem with a nonlinear dynamical 
constraint as an equivalent problem in which the system model was linear. The 
solution to the linear system was approximated using collocation techniques, and 
the optimal control problem approximated as an unconstrained nonlinear program. 
In particular, the optimal momentum management problem was reduced to an 
unconstrained, nonlinear, quadratic programming problem, and a computational 
algorithm for obtaining solutions was described. 

A solution to the optimal momentum management problem yields a state and 
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control history that maintains both the spacecraft attitude and momentum near 
desired states in the presence of external disturbances. The primary advantage of 
the proposed approach is that issues relating to attitude control and momentum 
management are treated in a combined framework that incorporates a nonlinear 
vehicle model. The importance of treating the nonlinear vehicle dynamics was 
highlighted in a number of examples. 

The first two examples using the optimal momentum manager showed how the 
solutions manipulated the total momentum of the vehicle/CMG system to reduce 
the influence of external torques and Euler coupling on the vehicle. A third exam- 
ple included an external torque roughly representative of those the proposed NASA 
space station may encounter on-orbit. The last two examples illustrated how opti- 
mal solutions could exploit the gravity gradient torque to eliminate net momentum 
stored in the CMG system. In all but the first example, the variation of vehicle 
momentum was of the same order-of-magnitude as the nominal values. These vari- 
ations are large enough to invalidate the small state excursion assumption made in 
the developement of existing schemes for momentum management. 

This thesis has investigated the application of feedback linearization theory to 
two problems arising in the control of spacecraft with CMGs. While the approach 
to CMG steering developed in Chapter 3 is no more computationally intensive 
than existing algorithms, the computational overhead associated with the optimal 
momentum management scheme presented in Chapter 4 can be significant. This 
burden is offset by the fact that the computational algorithm developed here is very 
robust, thus fairly reliable estimates on convergence times can be made. The reliable 
convergence features of this approach to momentum managment make the technique 
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suitable for autonomous implementation when turn-around time is not a critical 
issue. In the space station examples studied here, the momentum management 
algorithm must be scheduled for execution about once an orbit (i.e. 90 min) 

Chapter 4 concluded with a description of a proposed spacecraft attitude control 
system based on the feedback linearized dynamics of the spacecraft. Although the 
detailed design of the closed-loop control system has been left for future research, 
the roles of the optimal momentum manager and the approximation algorithm de- 
veloped in Chapter 3 were emphasized. The optimal momentum manager played 
the role of a command generator, and would typically be called once per orbit, or 
when needed. The approximation algorithm was used to distribute environmental 
disturbances among the available CMGs. 

One limitation inherent in the proposed approach to momentum management 
is that it does not allow for discontinuous disturbances or control action. Step or 
impulse disturbances represent an important class of torques associated with dock- 
ing or configuration changes and are important to consider. Although high-order 
approximations can be used to deal with such disturbances, this is an inefficient 
approach. One possible strategy might be to use polynomial approximations on in- 
tervals that have their endpoints at the times when the discontinuities occur. Note 
that this technique could also be used to include actuators like jets in the momen- 
tum management problem, since their net effect on momentum appears linearly in 
both the nonlinear and the linearized models. 

Future investigations might also consider methods of establishing an optimal 
trajectory that would be valid on succesive orbits. One approach would be to apply 
the present algorithm to a time interval containing several orbits, and choose a 
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suitable subinterval of this solution. Another possibility would be to include the 
initial and final states as decision variables in the optimization, and to require 
that certain compatability conditions between the initial and final states be met. 
Such a technique should also eliminate the errors resulting from the polynomial 
approximation that were seen to occur near the final time in the examples that 
have been presented. 

In the initial stages of this research, it was found that a criterion that mea- 
sured the optimality of a CMG configuration could be chosen as a state variable 
in the linearized model. A number of simulations were done in which the gain of 
the CMG configuration was a state variable of the linearized system. One difficulty 
with this approach was that the transformation would fail to exist at points where 
this function attained a local optimum. In addition, when the CMG system has 
a large number of degrees-of-freedom, it is necessary to create a number of cri- 
teria that must be controlled independently. Although there may be limitations, 
further investigation of this idea may lead to useful steering algorithms for CMG 
configurations with few redundant degrees-of-freedom. 

One important question left for future investigation involves how uncertainty in 
the vehicle model effects the stability of closed-loop control systems designed using 
feedback linearizing transformations. Results are available which guarantee that 
sufficiently small errors will not destabilize such control systems [56]. For practical 
purposes, it would be useful to have quantitative bounds on the uncertainty that 
these systems can accomodate. Another related issue is that the linearizing trans- 
formations assume perfect state information. In practice, all state variables may 
not be available, and this may pose an additional limitation in certain applications. 
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The manner in which the linearizing transformations assimilate perfect model 
information might also be exploited in an adaptive control scheme. In applications 
such as the proposed NASA space station, where the spacecraft will be erected over 
a long period of time, the performance robustness of any fixed control system is a 
delicate question. If accurate model information can be provided to the algorithm 
constructing the linearizing transformation, the transformations will adapt perfectly 
to the changing system. 

In conclusion, this thesis has used feedback linearization techniques to develop 
approaches to the resolution of redundancy in actuator systems, and the solution 
of optimal control problems. Both of the problems studied in this thesis arise 
in various forms in other areas of control and systems theory. In particular, the 
similarity of spacecraft control with CMGs to the control of robotic manipulators 
has been discussed in some detail. The ideas developed in this thesis could similarly 
be applied to obtain practical solutions to problems in a number of related fields. 
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Appendix A 
Notation 


h 


i 

In 

M,N 

dim(-) 


q 


q* A 


qi 


qL 

LVLH 

R n 

X) /) • • • 


X 


denotes a momentum vector 
denotes an inertia tensor 

denotes an identity matrix of dimension n x n 
used to denote a manifold 
dimension of a linear vector space 
a quaternion 

quaternion of reference frame B wrt. A 

quaternion of a body fixed frame wrt. an inertial frame 

quaternion of a body fixed frame wrt. an LVLH frame 

Local- Vertical/Local Horizontal 

the real n-tuple space 

lowercase italics are the workhorses meaning is defined 
by context, e.g. n-tuple, function, etc. 
time derivative of x 
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A-b 


I */ 

x 

(U,t) 

A 

9 

U 

O 

T 

min ./(a) 
a 

max J(a) 
a 


dot used to emphasise an arithmetic product or matrix multiplication 

unit vector, e.g. h 

denotes the dual of a vector space 

denotes the Lie derivative with respect to the vector field / 

denotes the annihilator of a distribution 

a coordinate system with coordinate chart <f) defined on U 

defined equal to or defined by 

Kronecker delta, 6^ = 1 if i = j, Sij = 0 if i ^ j 

calligraphic letters denote distributions 

used to denote a set of CMG gimbal angles 

denotes and angular rate 

composition, / o g = f(g) 

denotes transpose 

minimize J{a) with respect to a 

maximize J(a) with respect to a 



Lie bracket of the vector fields / and g 
n-fold Lie bracket ad^g = g , ad^g = [/, a<Ff~ l g] 
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Appendix B 

Solution of the Feedback Linearization Problem 


In this appendix the external linearization problem is precisely defined and solved. 
Before discussing the problem, requisite background concepts and definitions are 
presented. The problem solved here is the same as that solved in [19]. The derivation 
presented here, however, follows the geometric approach of [57]. Conditions are 
stated which guarantee the existence of coordinate changes in the state and input 
spaces which together with nonlinear feedback linearize the state equations of a 
locally controllable nonlinear system. 

B.l Background from Differential Geometry 

In this section, the notation and concepts used to derive the results of subsequent 
sections are presented. In the next subsection, the notion of a differentiable manifold 
is introduced, along with a number of structures defined on manifolds. A completely 
integrable manifold is defined and the Frobenius Theorem is stated and interpreted. 
Most of the definitions and results can be found in greater detail in texts such as 
[58], [59], and [60], or in the appendix of [57]. 

B.1.1 Manifolds, Vector Fields, and Covector Fields 

The notation / : A — » R n denotes a vector valued function defined on an open set 
A C R m taking values in the real n-tuple space. Such a function is called smooth, 
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differentiable, or of class C°° if all partial derivatives of all orders of the component 
functions /i,. . . ,/ n exist and are continuous. Nearly all of the results which follow 
are true assuming that all n th -order partial derivatives of the component functions 
are continuous for n sufficiently large. The derivative (or Jacobian) of / is denoted 
D/(*) and is given by the n by m matrix of partial derivatives 


P>/(*)]y 


dft 


[dxj\ 


(*) 


(B.l) 


In the case where m = n, and the matrix D/(®) is nonsingular at x, f is called 
a diffeomorphism. By the Inverse Function Theorem, differentiability of / and 
nonsingularity of (B.l) at x guarantee that the function / has a smooth inverse on 
an open neighborhood of f{x). 

A smooth k-manifold M in R n is a set of points in R n having the property: 
Given p € M, there are open sets A, B in R n containing p, and a diffeomorphism 
g : A — » B such that 


g [A n (R fc x 0"- fe )) = B n M (B.2) 

The notation X x Y denotes the cross product of the sets X and Y ; this is the set 
of ordered pairs of the form [x,y] where x £ X and y £ Y. In particular, the set 
R* x 0 n ~ fc denotes the set of n-tuples for which the last n — k elements are zero. 
This definition of a manifold is more restrictive than usually employed in abstract 
discussions; the structure just defined is often called a regular submanifold of R n . 

From a practical point of view, most dynamical systems can be thought of 
as evolving in an n-dimensional state space , and since the linearization theory to 
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Figure B.l: A Two-dimensional Manifold in R 3 

be developed is local, it is not necessary to have a definition which characterizes 
manifolds more generally. The definition of a manifold used in this discussion is 
illustrated in Figure B.l for a manifold consisting of a torus in three-dimensional 
space. The function defined by the rule 

a(x) = sea* (B.3) 

is called a coordinate patch on M, and its inverse <f)(p ) = a _ 1 (p) is called a coor- 
dinate chart on M. Note that because the manifold is a subset of R n , the range 
of the coordinate patch a is an n-tuple with n > fc. Since the derivative Da is 
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not a square matrix, the mapping a will not generally have an inverse on R n . The 
coordinate patch does have an inverse on an open set of points in the manifold. 
Similar remarks apply to the coordinate chart which is a map <f> : R n — > R fe . 

Based on these definitions, a manifold can be envisioned as a smooth surface 
in R n , and a coordinate patch is an invertable function that maps points in a 
^-dimensional Euclidean space onto the manifold, M. The coordinate patch and 
coordinate chart can be thought of as defining the relationship between points on 
a flat map covered with a rectangular grid and a spherical globe. The definition 
says that for each point on the manifold, p, there is an open set U C M containing 
p, and a coordinate chart (f> defined on U; the pair ( U,<f > ) is called a coordinate 
system at p. The principal advantage of defining a manifold in a Euclidean space is 
that functions defined on the manifold can be extended to functions defined on the 
surrounding space in a natural way. This will simplify notation when it is necessary 
to manipulate functions defined on the manifold; definitions can be expressed in 
the local coordinates of the Euclidean space surrounding the manifold rather than 
in the coordinates of the range space of the coordinate chart. 

A tangent vector v , with components {vi, . . . ,v n } at a point x in R n , can be 
thought of as an operator on differentiable functions. If A : R n — > R is a smooth 
function on a neighborhood of x, then u(A) is the real number which is the directional 
derivative of A in the direction of the vector {i^, . . . , v n }. In the standard coordinate 
system for R n - which is itself a manifold - v(A) is given by 

«»-!>(£). (B - 4) 
Denoted T x R n , the tangent space to the manifold R n at x is the set of all tangent 
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vectors at x. The tangent space is a real, linear vector space with scalar multipli 
cation and addition defined by 


(r-.ii, = r, 




The set of tangent vectors which are directional derivatives along each of the n 
coordinate axes in R n are denoted by 



(B.6) 


and form a basis for the space T x R n in much the same way that the vectors e; = 
[Oi, . . . , lj, . . . , 0 n ] form a basis for R n . In fact, the e* are the components of the 
basis vectors of the tangent space; consequently, the vector space stucture of T x R n 
can be envisioned with the same geometric structure as the vector space R n . It 
is possible to define the tangent space at a point on a fc-manifold by introducing 
a coordinate system, (C/ ? 0), at a point p £ M. Denoted T P M, this is also a real 
vector space. Associated with the coordinate chart, there are k tangent vectors at 
p defined by 



A 


A 



[DA][D^r‘].=*(p) 


(B.7) 


where D<^ 1 denotes the I th column of the derivative of <f> 1 and y denotes a point in 
the space containing the manifold. If <f>~ 1 is imagined as mapping the i ih coordinate 
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axis onto a parameterized curve on the manifold, the i th basis vector takes the 
derivative of A in the direction of the tangent to the curve at p. In other words, 
if <£ -1 is viewed as a function of a single variable, say t = x,-, then the tangent 
to the curve can be interpretted as a velocity vector. In local coordinates, this 
velocity vector is given by D^f 1 . Since <f> is a diffeomorphism of open sets, the 
set of directions associated with the k tangent vectors (B.7) are independent and 
lie in a fc-dimensional planar surface tangent to M at p. These k tangent vectors 
form a basis for the tangent space of M at p called the natural basis induced by the 
coordinate system. Every element of T P M can be expressed in the form 


* - i -■(&), (b - 8) 

If A is a smooth function defined on an open neighborhood of p G M, and v G T P M , 
then v(A) is the derivative of A in a direction which is tangent to the surface of M 
at p. For manifolds in R n , it is useful to associate a tangent vector with a direction 
in R n and a point on the manifold. 

Let (U,<f>) and (V^) be coordinate systems about a point p on an m-manifold 
M. The two coordinate charts induce different bases for the same tangent space 
T P M. The relationship between the two bases can be found by expressing each basis 
vector (^r)p(A) in terms of the basis induced by (U } <j>). 

( d(X o (j)~ l o <f> o 

\ / y=x p(p) 

a / g(Aon ] 

y=x V d x j J I= ^( p ) V dyi ) y =v.(p) 
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E 

;= i 


' Who*- 1 ) ' 
k <9y; 


y=^(p) 



(B.9) 


The function defined by / = 0 o i/> _1 is a diffeomorphism and defines a change of 
coordinates about p £ M. If {u x , . . . , u m } and {ttfi, . . . ,u; m } denote the components 
of the same tangent vector in the bases induced by (U,<j>) and (y, i/>), respectively, 
then (B.9) and the requirement that 




1 = 1 


Mi, 


3=1 




(B.10) 


imply that 


> 

" 

Wl 


: 

► = 



J 
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d A 

dxi 

dfrn 

dxi 


dh_ 

dx m 

dfrn 

dx m 


" 

r 


Vi 

4 

• > 




[ = Df(x)v 


(B.ll) 


The notion of a change of basis can be generalized to an arbitrary mapping 
between two manifolds. Let N and M be smooth manifolds of dimension n and m, 
and let f : N — > M be a smooth mapping. The differential of / at p £ N is a 
linear map /* : T P N — > T/( P )M defined by 


(/.(«)) (*) = «(W) (B.12) 

for A is a smooth real-valued function. An analogous computation to (B.9) can be 
applied to v £ T P N and w £ T/( P )M to show that the component vectors are related 
by w = [D/(z)] v . 
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A vector field on a fc-manifold M is a map assigning a tangent vector to each 
point on M. A smooth vector field v(q) is such that for each point p, there is a 
coordinate chart ( U , <j>) about p, and k real-valued smooth functions vi(q), . . . ,vi,(q) 
defined on U such that 

v(i) = E »<(«) (4;) ? < B - 13 ) 

for all q (E U. It is often useful to think of the vector field as defined by the set 
of component functions {^( 9 )} and to associate the vector field with the column 
vector of component functions taking the form 


v(q) = 


v i(?) 


V M 


(B.14) 


It is a customary abuse of terminology to refer to the vector of component functions 
as a vector field; it is more accurate to say that they define a vector field. If the 
manifold under consideration is R n , the vector of n component functions can be 
naturally associated with a direction in the space. 

The Lie derivative of a scalar valued function A with respect to a vector field is 
denoted L V A and defined to be 


(L„A)(p) = (»(p))(A) 


(B.15) 



The Lie derivative can be expressed in local coordinates yielding 


(L,A)(p) 


dx ax 

dxi dx r 




(B.16) 


If one of the component functions of the coordinate chart is used in place of A, then 
(B.16) implies that L v <j>i is the projection of the components of v on the gradient 
of </>j. This implies that (B.13) can be written 


»(,) = £ L,«,) (^-) (B.17) 

The vector space of all smooth tangent vectors on a manifold can be given the 
structure of a Lie algebra by defining the Lie bracket of a pair of vector fields. The 
Lie bracket of two vector fields v, w is another vector field denoted [u, tn] and defined 

by 

(M(P))(A) = L„(L„(A))(p)-ML.(A))(p) (B.18) 


In local coordinates, direct computation shows that the components of the vector 
field [v,w] are given by 


[v,u;](z) = [Dtu(x)]{n(a3)} — [Dn(®)]{u?(a:)} (B.19) 

It is instructive at this point to examine the connection between vector fields, 
as defined here, and differential equations. It is customary to consider a system of 
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n autonomous equations of the form 


^jT = /i(*i(0i- ••>*«(<)) 

: (B.20) 

= /n(«l(<)r",®nW) 

A parameterized curve x(t) € R n is a called a solution if it satisfies (B.20). Viewed 
as directions in R n , both sides of (B.20) can be associated with vector fields in a 
natural way, namely, 


/(*) = £/<(*. 


(B.21) 


and 




(B.22) 


Note that f(x) is being used to denote a tangent vector field in this context, the 
fi{x) are component functions. At any particular time, X(t) is also a tangent vector 
at the point z(i) to the 1-manifold which is the solution; the components ^ are 
merely the expression of X in the local coordinates of the manifold. Recalling the 
definition of the differential, (B.12), and regarding the solution x(t) as a mapping, 
x(t) : R 1 — > R n , it is possible to conclude that 


m = x. (| 


(B.23) 
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and since (B.21) and (B.23) must be equal along the solution, it is necessary that 


*.(<) (!) = /(*«) (B.24) 

This leads to the general notion of an integral curve of a differential equation on a 
manifold, M. A smooth curve, <r : (to,ii) — * M is an integral curve of the vector 
field / defined on the M if 

T.(|) = /MO) (B.25) 

for all t in the open interval (t 0 , <i). Clearly, an integral curve is simply a generalized 
notion of a solution to (B.20). The advantage of this description is that it is a 
coordinate-free notion of a differential equation (i.e. state variables have neither 
been defined nor chosen). 

The dual space to T p M is called the cotangent space and is denoted T*M , the 
tangent covectors are linear operators on tangent vectors that take values in the 
real numbers. The value of v * £ T*M at v £ T P M is denoted (v*,v). The space 
T*M has the structure of a real vector space with addition and scalar multiplication 
defined by 

(ciU* +€!«*,«) = Ci{T;*,u) + c 2 (v* 2 ,v) , C!,C 2 eR (B.26) 

For a smooth function A : M — ► R, there is a way of associating the differential 
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A* at p with a covector. Since A*(v) can be written in the form 


A.(v) = c(v)^ 


(B.27) 


and since A* is a linear operator, the scalar function c(v) must be a linear function 
of v. Since c(v) is a linear operator on tangent vectors, there must be a covector 
(dX) p such that 

A.(») = ((")„”)£ (B. 28 ) 


In local coordinates, the component in this expression is 


((d\) p ,v) = 


d\ 

dxi 


8X 

dxi 


*1 


(B.29) 


The last equality implies that 


{d<t>\ 


(~ 


d 

dh 


(*) 


dxj 

dx i 

dxj 


(B.30) 


r 

0 if i ^ j 

1 if i = j 
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and the covectors {d<f > 1} . . . , d<f> n } form a basis for the cotangent space which is dual 
to the basis for the tangent space. Any cotangent vector can be represented in this 
basis by an expression of the form 

«• = £>?(<%), < B - 3i > 
t=l 

Note that there are natural bases for the tangent space and the cotangent space 
induced by a coordinate chart on the manifold. If v and v* are expressed in the 
natural basis induced by the same coordinate chart at p, then 

( v % v ) = Jl v i{p) ■ v i(p) ( B - 32 ) 

i=i 

and it is natural to associate the operation of a cotangent vector on a tangent vector 
with an inner product on R fc . 

As in the case of a vector field, is is possible to define a covector field on a 
fc-manifold. A covector field on a manifold is a map assigning to each point of the 
manifold a cotangent vector. A smooth cotangent field u*(q) is such that for each 
point p G M, there is a coordinate system (£/,<£) about p and k real- valued, smooth 
functions Ux(g), • • * ^(g) defined on U such that 

v *(g) = ]£ v i(q)( d <f>i)p (b.33) 

i=1 

for all q G U. Smooth covector fields are also known as 1-forms. The notation d 
is a standard notation for the differential operator on fc-forms. Since a real- valued 
function is often defined to be a 0-form, (B.29) indicates that it is natural to think 
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of ( dX ) as a linear operator on tangent vectors; a 1-form. A smooth covector field 
is called an exact one- form if v*(p) = dX for some smooth real valued X. Later in 
this chapter, it will be shown that the linearization problem is solved by a certain 
set of exact 1-forms; the computational problem is to find the associated 0-forms. 
Allowing for a slight abuse of terminology, a row vector of n real functions will often 
be referred to as a covector field. If the manifold happens to be R n , the operation 
(v*,v) can be thought of as an inner-product of vector valued functions, u*(x) v(x). 

In proving the theorems in the next section, there will be a need for a Leibnitz- 
type formula relating Lie derivatives of functions, Lie brackets (which are Lie deriva- 
tives of vector fields), and Lie derivatives of 1-forms. Let v*(x) denote the (row) 
vector of component functions of a covector field v *, and let u(x) denote the (col- 
umn) vector of component functions of a vector field v. The Lie derivative of v* with 
respect to v is a covector field denoted L„v*. In local coordinates, the components 
of L V v* are given by the row vector of scalar functions 

L„u*(x) = u(x) t [Di;*(x)] + v*(x)[Du(x)] (B.34) 

Let / and g denote smooth vector fields and let v* denote a smooth covector field, 
the Leibnitz-type formula that will be needed is 

L /(u*,g) = (LfV*,g) + (v*,[f,g]) (B.35) 

If v* is an exact one-form, then v* = dX for some real- valued function A, and (B.35) 
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can be rewritten in the two equivalent forms 


L ,(d\,g) = (L,(iX),g) + («,(/, j]) (B.36) 

L,l,\ = « X/A),») + <<fA, [/,«,]> (B.37) 

The two equalities (dA,/) = L/A and L/dA = d( L/A) are used to write (B.36) in 
the form (B.37). 

B.1'2 Distributions and the Frobenius Theorem 

For d < m, a smooth d-dimensional distribution V on an m-dimensional manifold M 
is a smooth mapping which assigns to each point p 6 M a subspace V{p) of T P M. A 
more insightful definition comes by considering a set of smooth, independent vector 
fields {ti , . . . , T^} defined on M. These vector fields define a distribution on M 
according to the rule 

V = span{ri : iG[l,...,d]} (B.38) 

The notation span{...} denotes the subspace of a vector space which is the set of 
all linear combinations of elements of the argument list. Just as a tangent vector 
can be associated with a vector in a Euclidean space, V(p) can be associated with 
a subspace of Euclidean space with origin fixed at p. The subspace spanned by the 
set of directions {r;} is the set of directional derivatives associated with elements 
in V(p). The dimension of a distribution V at p is the dimension of the subspace 
V(p) of T P M. A point p € M is a regular point of a distribution if there is an open 
neighborhood U C M containing p such that the dimension of the distribution is 
constant on U . A vector field is contained in a distribution on a manifold, written 
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r E V, if r(p) E V{jp) for all p E M. If V is a smooth d - dimensional distribution 
in a neighborhood U of a regular point p, then there is a set of d vector fields 
{t< : i E [l,...,d]} and d smooth functions {c,- : i € [l,...,d]} such that every 
r € T> can be represented in the form 

r = ^CjTi (B.39) 

i=i 

on U. The set of d vector fields on U is called a local generator for the distribution. 

A distribution is involutive if the Lie bracket of any pair of vector fields is 
contained in the distribution. The statement that a distribution is closed under the 
bracketting operation is equivalent to requiring that the Lie bracket of two vector 
fields contained in the distribution can be expressed locally in the form of (B.39). If 
/ is a vector field and V a distribution on M, then [/, V) is the distribution defined 
by the rule 

[ /.®] = {[/.r]: Vr 6 0} (B.40) 

It is natural to also introduce the objects which are dual to distributions. A 
codistribution T>* on a manifold M is a mapping assigning to each point of M a 
subspace of T*M . If V* is a smooth codistribution of dimension d at p, then there 
is a set of d independent covector fields r* defined on a neighborhood U containing 
p such that 

V* = 3pan{r* : (B.41) 

Associated with a distribution D is a codistribution called the annihilator of D, 
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denoted V s - , and defined by the rule 


V x = {v* € T;M : (v» = 0 for all v G V(p)} (B.42) 

If V is a d- dimensional distribution on an n-dimensional manifold, then dim(V) + 
dimiV 1 ) = n. In the Euclidean domain, this statement says that the dimension of 
a subspace plus the dimension of the orthogonal subspace is equal to the dimension 
of the entire space. 

A nonsingular d-dimensional distribution V on an n-manifold N is completely 
integrable if at each p G N there exists a coordinate chart ( U,<j > ) with coordinate 
functions <^i, . . . , such that 

%)=^n{(A) f ,...,(^)J (B.43) 

for all q E U. If there exists a coordinate chart such that (B.43) is satisfied, 

then it is possible to partition U into slices passing through any p 6 U defined by 

S p = {q G U : <f*i(q) = for d + 1 < i < n} (B.44) 

The slices S p are d-dimensional manifolds which have the property that their tangent 
spaces coincide with the subspaces V(p) at each point of U. Another important 
aspect of completely integrable distributions is that they are involutive. This can 
be seen by considering the vector of component functions of an arbitrary tangent 
vector t(;c) = {r x (aj), . . . ,r n (x)} T in a completely integrable distribution. Definition 
(B.43) implies that there is a change of coordinates under which the last n — d 
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component functions are identically zero; therefore, the last n — d rows of Dr(x) 
are identically zero. Applying (B.19) to two such vector fields implies that the last 
n — d component functions of the bracket are also zero and the bracket is contained 
in the distribution. 

Clearly involutivity is a necessary condition for a distribution to be completely 
integrable; it is in fact true that it is also a sufficient condition. 

Frobenius Theorems A nonsingular distribution is completely integrable if 

and only if it is involutive. 

The proof of sufficiency is fairly involved and can be found in any of the references 
cited in the introduction to this section. 

The Frobenius Theorem leads to an algorithm for the construction of a change 
of coordinates that satisfies (B.43). Let V be a nonsingular involutive d-dimensional 
distribution, a coordinate system about p £ N, and {r; : i € 1, . . . , d} a local 

generator for V on U. In the natural basis of the chart, the r,- can be represented 
by 

Ti = (B ' 45) 

The objective is to find a new coordinate system ( V about p such that V == 
span{ . . . , In this basis, the set of local generators is given by 

T. = E(i^)( a -~) (B ' 46) 
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Since the coordinate transformation must satisfy (B.43), it is necessary that 


L = 0 for 1 < i < d and d + 1 < j < n 


(B.47) 


In order to insure that the manifold have the same rank in the new coordinates, the 
matrix defined by 

r 

I'nV’i ... L Tl V>d 


K = 


(B.48) 


I-'TiV’l • • • 

must have rank d on V. The coefficients of the d generators on V are given by 


M x ) = it forrfk) O <f> *(x)] d ^ - 3 - ° d x ' k ~ 


k= 1 

n 




fc = 1 


dil 

dx k 


= (D«*)][T(*)] 


(B.49) 


where £ has been introduced to denote the change of coordinates on U. The last 
equation, together with the requirements (B.47) and (B.48), yield the system of 
linear partial differential equations 


[D«x)] [T(x)] 


K(x) 

n-d)xd 


(B.50) 


whose solution yields the change of coordinates function. 

A solution of (B.50) has been interpreted as defining a transfomation under 
which a distribution has certain properties; there is a dual interpretation of the 
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equation which indicates that solutions also characterize a certain codistribution. 
Solving the system of partial differential equations (B.50) is equivalent to finding 
rt — d real-valued functions . ..,(„ona neighborhood V of p with the property 
that the covectors d^ +1 (p), . . . , d£ n (p) are linearly independent and 

{d£i {p)>Tj(p)) — 0; for d + 1 < i < n and 1 < j < d (B.51) 

for all p £ V. The conditions (B.51) imply that the set of n — d covector fields {d&} 
are linearly independent on V , and span the (n — d)-dimensional codistribution V L . 
This interpretation of the Frobenius theorem will be used extensively in the proof 
of the linearization theorem stated in the next section. 


B.2 Solution of the Feedback Linearization Problem 

In this section, the linearization problem is solved for the nonlinear system of equa- 
tions 

m A 

£ = /(*) + 22ffi( x ) u i — f( x ) + G(x)u (B.52) 

1=1 

where /, . . . ,£r m are smooth vector fields, x £ R n , and u £ R m . The goal is 

to construct a state transformation z = F(x) and an input transformation u = 
a(x) + /3(x)v such that the system associated with the pair (z,v) in Figure 2.4 
is linear. Two systems related in this way will be called feedback equivalent . For 
the transformation to be invertible, it is also necessary that the matrices /3(x) and 
D F(x) be nonsingular. Upon substitution into (B.52), the transformations (jF, a,/?) 
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must satisfy 


i = [DFo(F-'( 2 »] 

m 

= Az + J2 M« 


(/ + 5Zs» a ») + 

»=1 t=l 


o (F-'(z)) 


i = 1 


= Az + Bv 


(B.53) 

(B.54) 

(B.55) 


In this section, the problem is solved assuming that {A,B) is a controllable pair. 

All controllable linear systems can be transformed using feedback and linear 
coordinate changes in the state and input spaces to a system in Brunovsky form 
[61]. This is a canonical form where the A matrix is block diagonal, and the i th 
block is a Kf x /c* matrix of the form 


0 1 0 ... 0 

0 0 1 ... 0 

0 0 0 ... 1 


0 0 0 ... 0 


(B.56) 
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The corresponding column of the 5-matrix is 



(B.57) 


where the nonzero entry corresponds the the last row of the submatrix A;. Note 
that all linear systems having a particular Brunovsky form are related by a trans- 
formation of the type being allowed for nonlinear systems, namely, feedback and 
coordinate changes in the state and input variable. If it is desired to solve the lin- 
earization problem for a particular pair, (A*f, 5^), it is sufficient to solve the problem 
for the pair (A, B) in Brunovsky canonical form. Without loss of generality, it can 
therefore be assumed that the system matrices (A, B) of (B.55) are in Brunovsky 
form. 

The Brunovsky canonical form is a multivariable generalization of the phase 
variable form for single-input linear systems, where feedback has been used to place 
all system eigenvalues at zero. Figure B.2 shows a realization of a multi- variable 
system in Brunovsky form, where the k* are assumed to be in decreasing order. The 
set of integers Ki correspond to the Kronecker indices of a minimal realization [37], 
Inspection of (B.53) and (B.54) indicates that the feedback linearization problem is 
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Figure B.2: Realization of a Brunovsky Canonical Form. 
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solvable if there exists a solution to the system of partial differential equations 


[D-F) • [/ + 2>“il 


t = l 


o(^(z)) 


[DF] • [G/3]i o (F -1 (z)) 


Az 

b{ for 1 < i < m 


(B.58) 

(B.59) 


on an open neighborhood about some point. In this section, necessary and sufficient 
conditions axe derived which guarantee the existence of solutions to these equations. 

Before proceding with the proof, it is necessary to state several definitions that 
will be used. Given the set of vector fields the set of distributions 

Qo,Gi,... are defined recursively by 


Go = 3pan{g l ,...,g m } (B.60) 

Gi = Gi- 1 + U,Gi- 1] (B.61) 

If each of the distributions Q{ is nonsingular, then 

Qi = 3pan{ad k fgj : 0 < fc < i, 1 < j < m} (B.62) 

where the operator ad k gj is the vector field constructed by the recursion adP^gj = gj 
and ad k gj = [/, adj^gj]. The importance of the result (B.62) is that it yields an 
algorithm for constructing a basis or local generator for the distribution It 
also yields an algorithm for constructing the dual basis for QP. In other words, 
if the vector fields ad k gj span an m-dimensional subspace of the tangent space to 
an n-manifold, there will be n — m exact covector fields {c/A>} - by the Frobenius 
Theorem - such that / ad k gj,d\i \ = 0 for each basis element of Qi. Since the one- 
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forms annihilate the basis elements of the distribution, every linear combination of 
the {dXi} annihilates every vector in the distribution Gi] therefore, the exact one- 
forms are a basis for Qt* The previous argument is used extensively in the proof of 
the linearization theorem. 

The sequence of integers Vo,Vi , . . . is defined by setting 

v 0 = dim Go (B.63) 

Vi = dim Qi — dim Gi-u i > 1 (B.64) 

(B.65) 

If denotes the last nonzero element in the sequence {vi}, and if dim Qi* = n, 
then 

i* 

y; Vj = n (B.66) 

1=0 

The set of indices {mo, . . . , m;-} is defined according to 

m o = v^ 
m 0 +mi = -i 

m 0 + mi + m 2 = v { *_ 2 

m 0 + h mi* = Vo 

If dim Gi m = 71 1 it follows that 

dim G^-i - i • mo + (i — 1) • m\ + 2 • m ;_ 2 + m;_i, 1 < i < i * (B.68) 


(B.67) 
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The last set of indices that is needed is given by 


Ki = t* + 1 

if 

1 < i < Vi* 

= i* 

if mi > 0 and 

Vi* + 1 < i < Vi*_ 

Ki = i* — 1 

if m 2 > 0 and 

Vi *- 1 + 1 < i < v^ 

K{ = 1 

if m;. > 0 and 

Vi + 1 < i < Vvg 


(B.69) 


The interpretation of the three sets of indices is made more clear by referring 
to Figure B.2. The number i* is the number of banks of integrators in the figure. 
The K{ are the number of integrators in the x th string of integrators, and hence, the 
Kronecker indices of the realization. The set of integers Vi denote the number of 
state variables that are integrated in the ( x + l) th bank of integrators - counting 
from the left. The m; are the number of state variables coming out of the (x + l) th 
bank of integrators - counting from the right - whose integrals are not state variables. 


The conditions under which the feedback linearization problem is solvable can 
now be stated [57]. 


Theorem: The feedback linearization problem is solvable on a neighborhood U 

of a point p if and only if p is a regular point of for all i > 0, 
dim(Qi*) = n, and the distribution is involutive for all x such that 
7^ 0- 


The requirement that p be a regular point of the distributions implies that the 
qualitative characteristics of the dynamical system are similar on a neighborhood 
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of p. If the differential equation (B.52) describes a physical system, the points which 
are not regular points often indicate pathological conditions (e.g. a change in the 
number of independent state variables). The requirement that dim{Gi • = n insures 
that the nonlinear system does have n independent state variables and is locally 
controllable. The involutiveness of certain Qi are the critical conditions which will 
be interpretted in the following proof. 

Without loss of generality, it can be assumed that the set of vector fields 
9u---i9m is independent, or equivalently, that v 0 = m. It will be seen from the 
construction which follows that the number of control variables necessary to solve 
the problem equals the dimension of the distribution Go- If one of the gi is not 
independent, it can be dropped by assuming U{ = 0 from (B.52) without affecting 
the construction of the linearizing transformation. The issue of redundant input 
variables will be discussed in more detail in Chapter 4. The sufficiency of the pre- 
ceeding theorem is proven in detail. Since the proof is constructive, the argument 
can be reversed to show necessity of the stated conditions. The proof of necessity is 
omitted since it is available in the literature, and since the construction used in the 
sufficiency proof is used in the next section to linearize the two systems of equations 
derived in the previous chapter. 

The most difficult part of the proof is the construction of a coordinate trans- 
formation about the point p. By assumption, the dimension of the distribution 
Qi • is n, and the distribution Gi '- 1 is involutive and has dimension n — m 0 ,' these 
assumptions imply that dim G^-i = rn<). From the Frobenius Theorem, there is a 
neighborhood U\ of p and m 0 exact one-forms on U x that span G^-i- Since the m 0 
covector fields are exact, each can be associated with one of m 0 smooth functions 
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/i 0 » such that 

(dh 0 i,v) = 0 for 1 < i < m 0 and V v € i (B.70) 

In particular, this is true for the elements ad^gj for 1 < j < m and 0 < k < i* — 1; 
consequently, the mo x m matrix 

M 0 (x) = [m?.(x)] = [(dhouatf gj) (x)] (B.71) 

has rank mo on JJ\. If this were false at some g, there would be a set of nonzero 
real numbers o 1? . . . , a mo for which 

X] °» ad l J gjj (g) = ( ^ o t - dh oi , ad'/gy Hg) = 0 (B.72) 

t=l \i = l / 

for all j £ [1 . . . m]. Taken with (B.70), (B.72) implies that there is a covector field 
c/A such that 

(dA,u) = ^aAiCg)-^?) = 0 VvG&.. (B.73) 

Since dim(^-) = n implies that dim((y t i) = 0, the covector dX = 0, and indepen- 
dence of the {dh 0 i} requires that a x = ••• = a mo = 0 contradicting the assumption. 

Now consider the distribution Gi *-2 of dimension n — 2m 0 — mi and the codis- 
tribution Gi ^-2 dimension 27710 + 771!. The Leibnitz formula (B.36) can be applied 
to the covector fields dh 0 i to show that 


- {dL } h 0U ad)g2j = (dh oi , ad k f +1 g^ - L f (dh oi , ad)g2j = 0 (B.74) 

for 1 < j < m, 1 < i < m 0 , and 0 < k < i* — 2 on the set U\. This necessarily 
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implies that 

{dLfdhoi, v) = 0 for 1 < i < m 0 and V v G £ t i_ 2 (B.75) 

It is also true that the set of covector fields 

dho\ 7 • • • ) d/iotnQ j dTj fh,Qi 7 . . . , dLy ho mo (B.76) 

are linearly independent for all x 6 U\. If the assertion were not true, there would 
be two sets of m 0 real numbers, {a t *} and {&{}, such that 

mo mo 

X) °* dh 0 i{q) + Y^ b idL f h oi (q) = 0 (B.77) 

1=1 1 = 1 

at some point q 6 This implies that 

0 = (]C£L°i a i dhoi(q) + bi dhfhoi(q) 7 ad l f l Qj S J 

{LTJx bi dhfhoi(q), a^^gj) (B.78) 

= - ( ET=\b i dh oi (q),a^g j ) 

for all 1 < j < m. The second equality follows from linearity of (•, •) and (B.70), 
and the third equality comes from the fact that hfdX = d(L/A) and the Leibnitz 
formula. The last equality can only be satisfied if all the bi are zero. Since bi = 0 
for all i and the d/i 0 » are independent, (B.77) can only be satisfied if a* = 0 for each 

i 

t. 

If m\ = 0, then we have constructed 2 m 0 vectors that span If m-i > 0, 

Qi * —2 Is involutive and the Frobenius theorem guarantees that there are rrii func- 
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tions h n (x), ... , hi mi such that the covector fields d/i n (i ), . . . , dh lmt , together with 
(B.76) are all linearly independent and 


( dhu,v ) = 0 for 1 <i < nil and V v € (?»•_ i (B.79) 


on some neighborhood [7 2 C U\ of p. 

The mi x m matrix Mi(x) is defined to be 

Mi(x) = [m,-j(x)] = ^dhu,ad l f ‘~ l gj^ (x)] (B.80) 

By construction, the (mo + mi) x m matrix given by 

Mo 

Mi 

has rank m 0 + mi at each point in the set 17 2 . If the statement were not true, there 
would be a set of real numbers such that at some point q £ {/ 2 

0 = (£2i a { dh 0 i, ad'j g ^ ( q ) + (E£L'i dh u , a<f f ~ l gj ) (q) 

= (- ETJ 1 + ££1 fc d/i H , (q) (B.82) 

= ^dA, ad'j~ l g^ (q) 

The second equality follows by the same trick used in (B.78). In combination with 
(B.82), the fact (B.79) implies that d\ must be contained in the annihilator of 
Since Qi*-\ is spanned by the covectors d/ion d\ must be a linear combination of 
them. Because of the independence of the covectors (B.76) and (B.79), (B.82) can 
only be true if all a; and b » are zero, a contradiction. 


(B.8i) 
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Continuing this procedure leads to a set of covector fields 

. . . , d/i>o mo > dljy ft-oi , . . . , dljy/iomo ? • • • ? dl/y hoi ? .*.jdIiy ^omo 
dh\\ , . . . , dh\mi ? • • • 5 dLy • • • > dliy hi mi 

• • • i 

There will be 

x*m 0 + (x* — l)m! +!••• + 2 t 7 V _ 2 + m i*-\ = n — m (B.84) 

of these covector fields, and by construction, the set is independent on some neigh- 
borhood U containing the point pJ Associated with each nonzero rrii there is an 
77i i x m matrix 

Mi = [m jfc ],(*) = [(dh ik ,ad t }~'g j S }] (®) (B.85) 

such that the matrix defined by 

(B. 86 ) 

has rank v x at all x 6 U . 

The covector fields given by 

dlij /loij • • - ? dliy hom 0 ? • • • 5 dLy , • • ■ , dL/y ^(t*— i)(mj* -1 ) (B.87) 
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together with (B.83) yield a set of (n — linearly independent one-forms. If 

there are in fact n such one-forms, the construction is complete; otherwise, m** > 0 
and there must exist a set of covector fields /ij-i, . . . , hi- m .. which combined with 
(B.83) and (B.87) span the n-dimensional cotangent space at p. 

At this point, the construction of a coordinate transformation is complete; the 
remaining step is merely to realize this fact. The n functions L have been 
constructed such that in local coordinates, their derivatives are independent about 
p. If these n functions are regarded together as a mapping on an n-dimensional 
manifold, the inverse function theorem guarantees that the mapping has a smooth 
inverse about p. This mapping is therefore a coordinate transformation about p. 
Consider one of the functions /iy(* ), what is the interpretation of the Lie derivative 
in the context of (B.52). Assuming that u = 0, L jhij(x) can be expressed 
in local coordinates as 


(B.88) 


. ^ dhij dx/g 

^ dx k dt 

Clearly, the Lie derivative Ly /i tJ (x) is the time derivative of /i,y(x) as the nonlinear 
system evolves without the influence of controls. More formally, the coordinate 
transformation is constructed by letting 


yn(x) = h oi (x) if 1 < i < 

yn(x) = hiu-v \(x) if mi > 0 and + 1 < i < i 

(B.89) 

Vii(x ) = if rrti • > 0 and u x + 1 < i < u 0 
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Figure B.3 shows the relationship between the functions constructed in the proof 
of the linearization theorem and the state-variables of the Brunovsky canonical 
form. The function yn(x) maps the state variable x to the output of the i th string 
of integrators in the Brunovsky form shown in Figure B.3. The state variables 
associated with the t th string of integrators of the system in Brunovsky form are 
given by 


»«(*) 

= h oi (x) 



ya(*) 

A 

= Vi i 

= L/h 0 i(*) 

(B.90) 

yui( x ) 

A . 

— 1) 

= L y-'M*) 


for 1 < i < m and the coordinate change in 

the state-space is 


2 = 


• ,ym(« m -l)(z)] T = F { x ) 

(B.91) 


The final part of the proof is to construct a linearizing feedback. Consider one 
of the coordinate transformations yn(x). From the construction of the coordinate 
transformation, the dual products 

(dL k f y iugj ) = 0 (B.92) 

for 1 < ijj < m and 0 < k < /c t - — 2. Together with the fact that (d\ 7 v) = L V A ? 
(B.92) implies that 

(dhjyn , gj'j = Jjg.h^yn (B.93) 

Because L ff L*yii(x) = 0 for 0 < fc < /c; — 2, the Lie derivative L*y*i(x) depends 
only on Jj^ +l yn{x). This is equivalent to L^L^y^x) = l; +i y tl (®) for the specified 
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Figure B.3: Canonical Form of Linearized System with State- variables. 
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range of k. The remaining step is to compute the feedback and the input trans- 
formation such that the control variables v* = LiLy i_1 y^(x) can be used to specify 
the derivatives of the variables y t ( H ._i). This is accomplished by first computing the 
value of V{ associated with a specified state-control pair (s, u) yielding 

IiiLy’ y,i(*) = I / /+Gul i y* y»i(®) 

= + iGuiy-'yaiv) g4 ^ 

= i»/*y<i(®) ! + i L giUi L K f '-'yii(x) 

= L K /ya(x) + EjLi (dL'^ <-1 yi 1 (a;),yj) u,- 

The m expressions for the individual components u; can be combined to obtain a 
single vector-matrix equation 


v 


t?!(. iM 


+ l ya(^)i9j}] u = ui(c) + j4(a;) u 


(B.95) 


By construction, the matrix A{x) is nonsingular on the neighborhood U containing 
p; therefore, (B.95) can be solved for the u which produces a desired v. Defining 


a(z) = —A x (x) 

u = a(x) + (3(x)v where 

(3(x) = A -1 (s) 

the feedback linearization problem is solved by the state coordinate change F(x), 
the state feedback a(x), and the linear change of coordinate in the input space (3(x ). 


L?yu(*) 


Ly" 


Vm l(x) 


(B.96) 
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